PWGHF/correlationHF/AliHFCorrelationFDsubtraction.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, 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 /* $Id: $ */
  17
  18 /////////////////////////////////////////////////////////////
  19 //
  20 // Class for subtracting D-hadron correlations from secondary D from B meson decay from D-hadron correlations of prompt+secondary D mesons
  21 //       n.b. requires the evaluation of the fraction of prompt D mesons using the same receipes used for D meson cross-section and RAA in D2H pag
  22 //       (-> fprompt obtained from FONLL predictions of prompt and secondary D, prompt and secondary D meson efficiencies and a range of hypotheses for
  23 //       RAA(DfromB)/RAA(promptD) (if needed)
  24 //
  25 // Author: A. Rossi, andrea.rossi@cern.ch
  26 /////////////////////////////////////////////////////////////
  27
  28 #include <Riostream.h>
  29 #include <TMath.h>
  30
  31 #include <TH1D.h>
  32 #include <TF1.h>
  33 #include <TFile.h>
  34
  35 #include <TH1D.h>
  36 #include <TGraphAsymmErrors.h>
  37 #include <TCanvas.h>
  38
  39 #include "AliHFCorrelationFDsubtraction.h"
  40
  41
  42 using std::cout;
  43 using std::endl;
  44
  45 ClassImp(AliHFCorrelationFDsubtraction)
  46
  47 //************** constructors
  48 AliHFCorrelationFDsubtraction::AliHFCorrelationFDsubtraction() : TNamed(),
  49   fptmin(),
  50   fptmax(),
  51   fhDataUncorrected(),
  52   fmethod(1),
  53   fgrConservativeFc(0x0),
  54   fgrConservativeNb(0x0),
  55   fnTemplates(0),
  56   fLastTemplAdded(0),
  57   fhTemplates(0x0),
  58   fhDataCorrected(0x0),
  59   fhRatioSameTemplate(0x0),
  60   fhRatioAllToCentralTempl(0x0),
  61   fcCompare(0x0),
  62   fcRatio(0x0),
  63   fCountTempl(0),
  64   fcAllRatio(0),
  65   fhEnvelopeMax(0x0),
  66   fhEnvelopeMin(0x0),
  67   fhEnvelopeMaxRatio(0x0),
  68   fhEnvelopeMinRatio(0x0),
  69   fgrEnvelope(0x0),
  70   fgrEnvelopeRatio(0x0),
  71   fSystUseRMSfromFlatDistr(1)
  72 { // DEFAULT CONSTRUCTOR
  73
  74 }
  75
  76
  77 AliHFCorrelationFDsubtraction::AliHFCorrelationFDsubtraction(const char* name, const char* title) :
  78   TNamed(name,title),
  79   fptmin(),
  80   fptmax(),
  81   fhDataUncorrected(),
  82   fmethod(1),
  83   fgrConservativeFc(0x0),
  84   fgrConservativeNb(0x0),
  85   fnTemplates(0),
  86   fLastTemplAdded(0),
  87   fhTemplates(0x0),
  88   fhDataCorrected(0x0),
  89   fhRatioSameTemplate(0x0),
  90   fhRatioAllToCentralTempl(0x0),
  91   fcCompare(0x0),
  92   fcRatio(0x0),
  93   fCountTempl(0),
  94   fcAllRatio(0),
  95   fhEnvelopeMax(0x0),
  96   fhEnvelopeMin(0x0),
  97   fhEnvelopeMaxRatio(0x0),
  98   fhEnvelopeMinRatio(0x0),
  99   fgrEnvelope(0x0),
 100   fgrEnvelopeRatio(0x0),
 101   fSystUseRMSfromFlatDistr(1)
 102 {// default constructor
 103
 104 }
 105
 106 AliHFCorrelationFDsubtraction::~AliHFCorrelationFDsubtraction(){
 107
 108
 109   delete fhDataUncorrected;
 110
 111
 112   delete fgrConservativeFc;
 113   delete fgrConservativeNb;
 114
 115   for(Int_t j=0;j<fnTemplates;j++){
 116     delete fhTemplates[j];
 117   }
 118   delete [] fhTemplates;
 119
 120   for(Int_t j=0;j<3*fnTemplates;j++){
 121     delete fhDataCorrected[j];
 122   }
 123   delete [] fhDataCorrected;
 124
 125   for(Int_t j=0;j<2*fnTemplates;j++){
 126     delete fhRatioSameTemplate[j];
 127   }
 128   delete [] fhRatioSameTemplate;
 129
 130
 131  for(Int_t j=0;j<3*fnTemplates-1;j++){
 132     delete fhRatioAllToCentralTempl[j];
 133   }
 134   delete [] fhRatioAllToCentralTempl;
 135
 136   for(Int_t j=0;j<fnTemplates;j++){
 137     delete fcCompare[j];
 138   }
 139   delete [] fcCompare;
 140
 141   for(Int_t j=0;j<fnTemplates;j++){
 142     delete fcRatio[j];
 143   }
 144   delete [] fcRatio;
 145
 146   delete fcAllRatio;
 147
 148
 149   delete fhEnvelopeMax;
 150   delete fhEnvelopeMin;
 151   delete fhEnvelopeMaxRatio;
 152   delete fhEnvelopeMinRatio;
 153   delete fgrEnvelope;
 154   delete fgrEnvelopeRatio;
 155 }
 156
 157
 158 Bool_t AliHFCorrelationFDsubtraction::Init(){// init histo, graphs and canvases
 159
 160   if(fnTemplates<=0){
 161     Printf("Number of templates not defined, aborting");
 162     return kFALSE;
 163   }
 164   fhTemplates=new TH1D*[fnTemplates];
 165   fhDataCorrected=new TH1D*[fnTemplates*3];
 166   fhRatioSameTemplate=new TH1D*[fnTemplates*2];
 167   fhRatioAllToCentralTempl=new TH1D*[fnTemplates*3-1];
 168   fcCompare=new TCanvas*[fnTemplates];
 169   fcRatio=new TCanvas*[fnTemplates];
 170
 171   fCountTempl=0;
 172   return kTRUE;
 173 }
 174 Bool_t AliHFCorrelationFDsubtraction::AddTemplateHisto(TH1D *h){
 175   if(fnTemplates<1){
 176     Printf("You first have to specify how many templates you want to use, aborting");
 177     return kFALSE;
 178   }
 179   if(fLastTemplAdded>=fnTemplates){
 180     Printf("You want to add more templates than those declared, aborting");
 181     return kFALSE;
 182   }
 183   fhTemplates[fLastTemplAdded]=(TH1D*)h->Clone();
 184   fLastTemplAdded++;
 185   return kTRUE;
 186 }
 187
 188 void AliHFCorrelationFDsubtraction::SubtractFeedDown(TH1D *hFDtempl){
 189   TGraphAsymmErrors *gr;
 190
 191   Double_t *xgr,*elxgr,*ehxgr,*fprompt,*elfprompt,*ehfprompt;
 192   if(fmethod==1||fmethod==3){
 193     gr=fgrConservativeFc;
 194   }
 195   if(fmethod==2){
 196     gr=fgrConservativeNb;
 197   }
 198   if(!gr){
 199     Printf("TGraph for feed-down subtraction not set. Aborting");
 200     return;
 201   }
 202
 203   xgr=gr->GetX();
 204   elxgr=gr->GetEXlow();
 205   ehxgr=gr->GetEXhigh();
 206   fprompt=gr->GetY();
 207   elfprompt=gr->GetEYlow();
 208   ehfprompt=gr->GetEYhigh();
 209
 210   Int_t binmin=-1,binmax=-1,bincenter=-1;
 211   Double_t ptcentr=0.5*(fptmin+fptmax);
 212   Double_t fpromptmin=1.,fpromptmax=0.,fpromptcentr;
 213   // Look for right bins
 214   for(Int_t j=0;j<gr->GetN();j++){
 215     if(binmin<0&&fptmin*1.00001>(xgr[j]-elxgr[j])&&fptmin<(xgr[j]+ehxgr[j]))binmin=j;
 216     if(binmax<0&&fptmax>(xgr[j]-elxgr[j])&&fptmax*0.9999<(xgr[j]+ehxgr[j]))binmax=j;
 217     if(ptcentr*1.00001>(xgr[j]-elxgr[j])&&ptcentr<(xgr[j]+ehxgr[j]))bincenter=j;
 218   }
 219
 220   if(TMath::Abs(xgr[binmin]-elxgr[binmin]-fptmin)>0.0001){
 221     Printf("Bin mismatch: xmin=%f; aborting ",xgr[binmin]-elxgr[binmin]);
 222     return;
 223   }
 224
 225   // central value for bin at centre of the bin (temporary)
 226   fpromptcentr=fprompt[bincenter];
 227
 228
 229   // Conservative approach: look for min and max values inside the pt range
 230   for(Int_t j=binmin;j<=binmax;j++){
 231     if(fprompt[j]+ehfprompt[j]>fpromptmax)fpromptmax=fprompt[j]+ehfprompt[j];
 232     if(fprompt[j]-elfprompt[j]<fpromptmin)fpromptmin=fprompt[j]-elfprompt[j];
 233   }
 234
 235
 236   TH1D *hFDCentr=(TH1D*)hFDtempl->Clone("hFDtemplCentr");
 237   hFDCentr->Scale(1.-fpromptcentr);
 238   printf("Centr fsec=%f \n",1.-fpromptcentr);
 239
 240   fhDataCorrected[fCountTempl*3]=(TH1D*)fhDataUncorrected->Clone(Form("hDataCorrectedTempl%dCentrFprompt",fCountTempl));
 241   fhDataCorrected[fCountTempl*3]->Add(hFDCentr,-1.);// Treatment of Stat unc.
 242   fhDataCorrected[fCountTempl*3]->Scale(1./fpromptcentr);
 243   fhDataCorrected[fCountTempl*3]->SetTitle(Form("Templ %d central f_{prompt}",fCountTempl));
 244
 245   TH1D *hFDmax=(TH1D*)hFDtempl->Clone("hFDtemplMax");
 246   hFDmax->Scale(1.-fpromptmin);
 247   Printf("Max fsec=%f ",1.-fpromptmin);
 248   fhDataCorrected[fCountTempl*3+1]=(TH1D*)fhDataUncorrected->Clone(Form("hDataCorrectedTempl%dMaxFprompt",fCountTempl));
 249   fhDataCorrected[fCountTempl*3+1]->Add(hFDmax,-1.);// Treatment of Stat unc.
 250   fhDataCorrected[fCountTempl*3+1]->Scale(1./fpromptmin);
 251   fhDataCorrected[fCountTempl*3+1]->SetTitle(Form("Templ %d min f_{prompt}",fCountTempl));
 252
 253
 254   TH1D *hFDmin=(TH1D*)hFDtempl->Clone("hFDtemplMin");
 255   hFDmin->Scale(1.-fpromptmax);
 256   Printf("Min fsec=%f ",1.-fpromptmax);
 257   fhDataCorrected[fCountTempl*3+2]=(TH1D*)fhDataUncorrected->Clone(Form("hDataCorrectedTempl%dMinFprompt",fCountTempl));
 258   fhDataCorrected[fCountTempl*3+2]->Add(hFDmin,-1.);// Treatment of Stat unc.
 259   fhDataCorrected[fCountTempl*3+2]->Scale(1./fpromptmax);
 260   fhDataCorrected[fCountTempl*3+2]->SetTitle(Form("Templ %d max f_{prompt}",fCountTempl));
 261
 262
 263   fcCompare[fCountTempl]=new TCanvas(Form("cCompareTempl%d",fCountTempl),Form("cCompare%d",fCountTempl),700,700);
 264   fcCompare[fCountTempl]->cd();
 265
 266   TH1D *hData=(TH1D*)fhDataUncorrected->DrawCopy();
 267   hData->SetTitle("Uncorrected");
 268   hData->SetYTitle("#frac{1}{N_{trig}}#frac{dN}{d#Delta#varphi} (rad^{-1})");
 269   hData->SetXTitle("#varphi_{ass}-#varphi_{trig} (rad)");
 270   hData->SetMarkerStyle(25);
 271   hData->SetMarkerColor(kBlack);
 272   hData->SetMarkerSize(1);
 273   hData->SetLineColor(kBlack);
 274
 275   fhDataCorrected[fCountTempl*3]->SetLineColor(kRed);
 276   fhDataCorrected[fCountTempl*3]->SetYTitle("#frac{1}{N_{trig}}#frac{dN}{d#Delta#varphi} (rad^{-1})");
 277   fhDataCorrected[fCountTempl*3]->SetXTitle("#varphi_{ass}-#varphi_{trig} (rad)");
 278   fhDataCorrected[fCountTempl*3]->SetLineWidth(2);
 279   fhDataCorrected[fCountTempl*3]->SetMarkerColor(kRed);
 280   fhDataCorrected[fCountTempl*3]->SetMarkerStyle(21);
 281   fhDataCorrected[fCountTempl*3]->SetMarkerSize(1);
 282   fhDataCorrected[fCountTempl*3]->Draw("sames");
 283
 284   fhDataCorrected[fCountTempl*3+1]->SetLineColor(kBlue);
 285   fhDataCorrected[fCountTempl*3+1]->SetMarkerColor(kBlue);
 286   fhDataCorrected[fCountTempl*3+1]->SetMarkerStyle(23);
 287   fhDataCorrected[fCountTempl*3+1]->SetMarkerSize(1);
 288   fhDataCorrected[fCountTempl*3+1]->Draw("sames");
 289
 290   fhDataCorrected[fCountTempl*3+2]->SetLineColor(kGreen-6);
 291   fhDataCorrected[fCountTempl*3+2]->SetMarkerColor(kGreen-6);
 292   fhDataCorrected[fCountTempl*3+2]->SetMarkerStyle(22);
 293   fhDataCorrected[fCountTempl*3+2]->SetMarkerSize(1);
 294   fhDataCorrected[fCountTempl*3+2]->Draw("sames");
 295
 296   fcRatio[fCountTempl]=new TCanvas(Form("fcRatioTempl%d",fCountTempl),Form("fcRatioTempl%d",fCountTempl),700,700);
 297   fcRatio[fCountTempl]->cd();
 298   fhRatioSameTemplate[fCountTempl*2]=(TH1D*)fhDataCorrected[fCountTempl*3+2]->Clone("hRatioMinFSecToCentr");
 299   fhRatioSameTemplate[fCountTempl*2]->Divide(fhDataCorrected[fCountTempl*3]);
 300   for(Int_t j=1;j<=fhRatioSameTemplate[fCountTempl*2]->GetNbinsX();j++){
 301     fhRatioSameTemplate[fCountTempl*2]->SetBinError(j,0.001);
 302   }
 303
 304   fhRatioSameTemplate[fCountTempl*2+1]=(TH1D*)fhDataCorrected[fCountTempl*3+1]->Clone("hRatioMaxFSecToCentr");
 305   fhRatioSameTemplate[fCountTempl*2+1]->Divide(fhDataCorrected[fCountTempl*3]);
 306   for(Int_t j=1;j<=fhRatioSameTemplate[fCountTempl*2+1]->GetNbinsX();j++){
 307     fhRatioSameTemplate[fCountTempl*2+1]->SetBinError(j,0.001);
 308   }
 309
 310   fhRatioSameTemplate[fCountTempl*2]->Draw();
 311   fhRatioSameTemplate[fCountTempl*2+1]->Draw("same");
 312
 313   fCountTempl++;
 314   delete hFDCentr;
 315   delete hFDmax;
 316   delete hFDmin;
 317
 318 }
 319
 320 TGraphAsymmErrors* AliHFCorrelationFDsubtraction::GetUncGraphFromHistos(TH1D *hRef,TH1D *hMin,TH1D *hMax){
 321
 322   //  Int_t npoints=hMin->GetNbinsX();
 323   Double_t ew=hMin->GetBinWidth(1)/2.;
 324   Double_t value,eyl,eym;
 325
 326   TGraphAsymmErrors *gr=new TGraphAsymmErrors();
 327   for(Int_t j=1;j<=hMin->GetNbinsX();j++){
 328     if(hRef){
 329       value=hRef->GetBinContent(j);
 330       eyl=hMin->GetBinContent(j)-value;
 331       if(eyl<0.)eyl*=-1.;
 332       if(hMax){
 333         eym=hMax->GetBinContent(j)-value;
 334         if(eym<0.)eym*=-1.;
 335       }
 336       else eym=eyl;
 337     }
 338     else {
 339       value=0.;
 340       eyl=hMin->GetBinContent(j)-1;
 341       if(eyl<0.)eyl*=-1.;
 342       if(hMax){
 343         eym=hMax->GetBinContent(j)-1;
 344         if(eym<0.)eym*=-1.;
 345       }
 346       else eym=eyl;
 347     }
 348
 349     gr->SetPoint(j-1,hMin->GetBinCenter(j),value);
 350     gr->SetPointError(j-1,ew,ew,eyl,eym);
 351   }
 352
 353   return gr;
 354 }
 355
 356
 357 TH1D* AliHFCorrelationFDsubtraction::ReflectHisto(TH1D *h,Double_t scale){
 358
 359   TH1D *h2=new TH1D(Form("%sReflected",h->GetName()),Form("%sReflected",h->GetName()),h->GetNbinsX()/2.,0.,TMath::Pi());
 360   for(Int_t j=1;j<=h->GetNbinsX();j++){
 361     Double_t x=h->GetBinCenter(j);
 362     Double_t y0=h->GetBinContent(j);
 363     Double_t ey0=h->GetBinError(j);
 364     Int_t j2;
 365     if(x>0&&x<TMath::Pi())j2=h2->FindBin(x);
 366     else if(x<0)j2=h2->FindBin(-1.*x);
 367     else if(x>TMath::Pi())j2=h2->FindBin(2.*TMath::Pi()-x);
 368     else {
 369       printf("Point %d excluded \n",j);
 370       continue;
 371     }
 372     Double_t y=h2->GetBinContent(j2);
 373     Double_t ey=h2->GetBinError(j2);
 374     h2->SetBinContent(j2,scale*(y+y0));
 375     h2->SetBinError(j2,scale*TMath::Sqrt(ey0*ey0+ey*ey));
 376   }
 377
 378   return h2;
 379 }
 380
 381 TH1D* AliHFCorrelationFDsubtraction::DuplicateHistoTo2piRange(TH1D *h,Double_t scale){
 382   if(h->GetBinLowEdge(h->GetNbinsX())+h->GetBinWidth(h->GetNbinsX())-h->GetBinLowEdge(1)>1.01*TMath::Pi()){
 383     Printf("AliHFCorrelationFDsubtraction: DuplicateHistoTo2PiRange works for histos with x range between 0 and pi");
 384     return 0x0;
 385   }
 386
 387   if(h->GetBinLowEdge(h->GetNbinsX())>1.01*TMath::Pi()){
 388     Printf("AliHFCorrelationFDsubtraction: DuplicateHistoTo2PiRange works for histos with x range between 0 and pi");
 389     return 0x0;
 390   }
 391   if(h->GetBinLowEdge(1)<-0.01){
 392     Printf("AliHFCorrelationFDsubtraction: DuplicateHistoTo2PiRange works for histos with x range between 0 and pi");
 393     return 0x0;
 394   }
 395
 396
 397   TH1D *h2=new TH1D(Form("%sTwoPiRange",h->GetName()),Form("%sTwoPiRange",h->GetName()),h->GetNbinsX()*2.,-TMath::Pi()/2.,1.5*TMath::Pi());
 398   for(Int_t j=1;j<=h2->GetNbinsX();j++){
 399     Double_t x=h2->GetBinCenter(j);
 400     Int_t j2;
 401     if(x>0&&x<TMath::Pi())j2=h->FindBin(x);
 402     else if(x<0)j2=h->FindBin(-1.*x);
 403     else if(x>TMath::Pi())j2=h->FindBin(2.*TMath::Pi()-x);
 404     else {
 405       printf("Point %d excluded \n",j);
 406       continue;
 407     }
 408     Double_t y=h->GetBinContent(j2);
 409     Double_t ey=h->GetBinError(j2);
 410     h2->SetBinContent(j,y*scale);
 411     h2->SetBinError(j,ey*scale);
 412   }
 413
 414   return h2;
 415 }
 416
 417
 418 TH1D* AliHFCorrelationFDsubtraction::GetHistoRelSystUncMin(){// Method to extrac the rel syst unc to be used later on in the analysis from the envelopes. IMPORTANT: remember
 419   // that with our procedure the rel. unc. is variation/distribution, and thus its value is strongly dependent on the baseline: only the absolute unc (->variation) is really meaningful
 420   // Options: if fSystUseRMSfromFlatDistr
 421   // =0 the full spread (->DeltaMax, DeltaMin) of the envelopes is used to define the syst unc,
 422   // =1 (default) DeltaMax/sqrt(3) and DeltaMin/sqrt(3) is used;
 423   // =2: as 1 but then the result is reflected to the range (0,pi), to decrease stat unc;
 424   // =3 as 1 then a smoothing is done (pol0 fit over 3 points), then the result is refelcted; other values not implemented yet (e.g. =2 could be taking the RMS)
 425   // =4 as 1 then a smoothing is done (linear fit over 3 points), then the result is refelcted; other values not implemented yet (e.g. =2 could be taking the RMS)
 426   if(!fhEnvelopeMinRatio){
 427     Printf("Cannot produce Rel syst unc histo: you first need to calculate the envelope");
 428     return 0x0;
 429   }
 430   TH1D *h=(TH1D*)fhEnvelopeMinRatio->Clone("hRelSystFeedDownMin");
 431
 432   if(fSystUseRMSfromFlatDistr==0){
 433     for(Int_t j=1;j<=h->GetNbinsX();j++){
 434       h->SetBinContent(j,fhEnvelopeMinRatio->GetBinContent(j)-1.);
 435     }
 436     return h;
 437   }
 438
 439   if(fSystUseRMSfromFlatDistr==1||fSystUseRMSfromFlatDistr==2||fSystUseRMSfromFlatDistr==3||fSystUseRMSfromFlatDistr==4){
 440     Double_t sqrtThree=TMath::Sqrt(3.);
 441     for(Int_t j=1;j<=h->GetNbinsX();j++){
 442       h->SetBinContent(j,fhEnvelopeMinRatio->GetBinContent(j)/sqrtThree-1./sqrtThree);
 443     }
 444     if(fSystUseRMSfromFlatDistr==1)return h;
 445
 446     if(fSystUseRMSfromFlatDistr==2){
 447       TH1D *hRefl=ReflectHisto(h);
 448       TH1D *hBackTo2Pi=DuplicateHistoTo2piRange(hRefl); // In this way we symmetrized to 0-pi
 449       delete h;
 450       delete hRefl;
 451       hBackTo2Pi->SetName("hRelSystFeedDownMin");
 452       return hBackTo2Pi;
 453     }
 454
 455     // case fSystUseRMSfromFlatDistr==3 or 4
 456     // Do smoothing and then symmetrize
 457     // For the smoothing, assume linear trend every 3 points, starting from transverse region
 458     Int_t jstart=h->FindBin(TMath::Pi()*0.5)-1;
 459     TF1 *f;
 460     if(fSystUseRMSfromFlatDistr==3){
 461       f=new TF1("mypol0","[0]",-TMath::Pi()*0.5,TMath::Pi()*1.5);
 462     }
 463     else if(fSystUseRMSfromFlatDistr==4){
 464       f=new TF1("mypol1","[0]+x*[1]",-TMath::Pi()*0.5,TMath::Pi()*1.5);
 465     }
 466
 467     TH1D *hOut=(TH1D*)h->Clone("hOut");
 468     for(Int_t j=jstart;j>=1;j--){// going towards 0
 469
 470       if(j>=3){
 471         Double_t xmax=h->GetBinLowEdge(j+1)+h->GetBinWidth(j+1);
 472         Double_t xmin=h->GetBinLowEdge(j-1);
 473         f->SetParameter(0,h->GetBinContent(j));
 474         if(fSystUseRMSfromFlatDistr==4)f->SetParameter(1,(h->GetBinContent(j)-h->GetBinContent(j-1))/(h->GetBinCenter(j+1)-h->GetBinCenter(j-1)));
 475         h->Fit(f,"REM","N",xmin,xmax);
 476         hOut->SetBinContent(j,f->Eval(hOut->GetBinCenter(j)));
 477       }
 478       else if(j==2){
 479         Double_t xmax=h->GetBinLowEdge(3)+h->GetBinWidth(3);
 480         Double_t xmin=h->GetBinLowEdge(1);
 481         Double_t y1=h->GetBinContent(1);
 482         Double_t y2=h->GetBinContent(2);
 483         Double_t y3=h->GetBinContent(3);
 484         h->SetBinContent(3,y2);
 485         h->SetBinContent(2,y1);
 486         h->SetBinContent(1,h->GetBinContent(h->GetNbinsX()));
 487         f->SetParameter(0,h->GetBinContent(2));
 488         if(fSystUseRMSfromFlatDistr==4)f->SetParameter(1,(h->GetBinContent(3)-h->GetBinContent(1))/(h->GetBinCenter(3)-h->GetBinCenter(1)));
 489         h->Fit(f,"RLEM","N",xmin,xmax);
 490         hOut->SetBinContent(j,f->Eval(h->GetBinCenter(2)));
 491         h->SetBinContent(3,y3);
 492         h->SetBinContent(2,y2);
 493         h->SetBinContent(1,y1);
 494       }
 495       else if(j==1){
 496         Double_t xmax=h->GetBinLowEdge(3)+h->GetBinWidth(3);
 497         Double_t xmin=h->GetBinLowEdge(1);
 498         Double_t y1=h->GetBinContent(1);
 499         Double_t y2=h->GetBinContent(2);
 500         Double_t y3=h->GetBinContent(3);
 501         h->SetBinContent(3,y1);
 502         h->SetBinContent(2,h->GetBinContent(h->GetNbinsX()));
 503         h->SetBinContent(1,h->GetBinContent(h->GetNbinsX()-1));
 504         f->SetParameter(0,h->GetBinContent(2));
 505         if(fSystUseRMSfromFlatDistr==4)f->SetParameter(1,(h->GetBinContent(3)-h->GetBinContent(1))/(h->GetBinCenter(3)-h->GetBinCenter(1)));
 506         h->Fit(f,"RLEM","N",xmin,xmax);
 507         hOut->SetBinContent(j,f->Eval(h->GetBinCenter(2)));
 508         h->SetBinContent(3,y3);
 509         h->SetBinContent(2,y2);
 510         h->SetBinContent(1,y1);
 511       }
 512
 513     }
 514     for(Int_t j=jstart+1;j<=hOut->GetNbinsX();j++){// now going towards 2pi
 515
 516       if(j<=hOut->GetNbinsX()-2){
 517         Double_t xmax=h->GetBinLowEdge(j+1)+h->GetBinWidth(j+1);
 518         Double_t xmin=h->GetBinLowEdge(j-1);
 519         f->SetParameter(0,h->GetBinContent(j));
 520         if(fSystUseRMSfromFlatDistr==4)f->SetParameter(1,(h->GetBinContent(j)-h->GetBinContent(j-1))/(h->GetBinCenter(j+1)-h->GetBinCenter(j-1)));
 521         h->Fit(f,"RLEM","N",xmin,xmax);
 522         hOut->SetBinContent(j,f->Eval(hOut->GetBinCenter(j)));
 523       }
 524       else if(j==hOut->GetNbinsX()-1){
 525         Double_t xmax=h->GetBinLowEdge(hOut->GetNbinsX())+h->GetBinWidth(hOut->GetNbinsX());
 526         Double_t xmin=h->GetBinLowEdge(hOut->GetNbinsX()-2);
 527         Double_t y1=h->GetBinContent(hOut->GetNbinsX()-2);
 528         Double_t y2=h->GetBinContent(hOut->GetNbinsX()-1);
 529         Double_t y3=h->GetBinContent(hOut->GetNbinsX());
 530         h->SetBinContent(hOut->GetNbinsX()-2,y2);
 531         h->SetBinContent(hOut->GetNbinsX()-1,y3);
 532         h->SetBinContent(hOut->GetNbinsX(),h->GetBinContent(1));
 533         f->SetParameter(0,h->GetBinContent(hOut->GetNbinsX()-1));
 534         if(fSystUseRMSfromFlatDistr==4)f->SetParameter(1,(h->GetBinContent(hOut->GetNbinsX())-h->GetBinContent(hOut->GetNbinsX()-2))/(h->GetBinCenter(hOut->GetNbinsX())-h->GetBinCenter(hOut->GetNbinsX()-2)));
 535         h->Fit(f,"RLEM","N",xmin,xmax);
 536         hOut->SetBinContent(j,f->Eval(h->GetBinCenter(hOut->GetNbinsX()-1)));
 537         h->SetBinContent(hOut->GetNbinsX()-2,y1);
 538         h->SetBinContent(hOut->GetNbinsX()-1,y2);
 539         h->SetBinContent(hOut->GetNbinsX(),y3);
 540       }
 541       else if(j==hOut->GetNbinsX()){
 542         Double_t xmax=h->GetBinLowEdge(hOut->GetNbinsX())+h->GetBinWidth(hOut->GetNbinsX());
 543         Double_t xmin=h->GetBinLowEdge(hOut->GetNbinsX()-2);
 544         Double_t y1=h->GetBinContent(hOut->GetNbinsX()-2);
 545         Double_t y2=h->GetBinContent(hOut->GetNbinsX()-1);
 546         Double_t y3=h->GetBinContent(hOut->GetNbinsX());
 547         h->SetBinContent(hOut->GetNbinsX()-2,y3);
 548         h->SetBinContent(hOut->GetNbinsX()-1,h->GetBinContent(1));
 549         h->SetBinContent(hOut->GetNbinsX(),h->GetBinContent(2));
 550         f->SetParameter(0,h->GetBinContent(hOut->GetNbinsX()-1));
 551         if(fSystUseRMSfromFlatDistr==4)f->SetParameter(1,(h->GetBinContent(hOut->GetNbinsX())-h->GetBinContent(hOut->GetNbinsX()-2))/(h->GetBinCenter(hOut->GetNbinsX())-h->GetBinCenter(hOut->GetNbinsX()-2)));
 552         h->Fit(f,"RLEM","N",xmin,xmax);
 553         hOut->SetBinContent(j,f->Eval(h->GetBinCenter(hOut->GetNbinsX()-1)));
 554         h->SetBinContent(hOut->GetNbinsX()-2,y1);
 555         h->SetBinContent(hOut->GetNbinsX()-1,y2);
 556         h->SetBinContent(hOut->GetNbinsX(),y3);
 557       }
 558
 559     }
 560     // Smoothing done, now simmetrize
 561     TH1D *hRefl=ReflectHisto(hOut);
 562     TH1D *hBackTo2Pi=DuplicateHistoTo2piRange(hRefl); // In this way we symmetrized to 0-pi
 563     delete h;
 564     delete hRefl;
 565     delete hOut;
 566     delete f;
 567     hBackTo2Pi->SetName("hRelSystFeedDownMin");
 568     return hBackTo2Pi;
 569
 570
 571   }
 572
 573   Printf("AliHFCorrelationFDsubtraction: wrong option for getting syst unc");
 574   return 0x0;
 575
 576
 577 }
 578
 579 TH1D* AliHFCorrelationFDsubtraction::GetHistoRelSystUncMax(){
 580   if(!fhEnvelopeMaxRatio){
 581     Printf("Cannot produce Rel syst unc histo: you first need to calculate the envelope");
 582     return 0x0;
 583   }
 584   TH1D *h=(TH1D*)fhEnvelopeMaxRatio->Clone("hRelSystFeedDownMax");
 585   if(fSystUseRMSfromFlatDistr==0){
 586     for(Int_t j=1;j<=h->GetNbinsX();j++){
 587       h->SetBinContent(j,fhEnvelopeMaxRatio->GetBinContent(j)-1);
 588     }
 589     return h;
 590   }
 591
 592   if(fSystUseRMSfromFlatDistr==1||fSystUseRMSfromFlatDistr==2||fSystUseRMSfromFlatDistr==3||fSystUseRMSfromFlatDistr==4){
 593     Double_t sqrtThree=TMath::Sqrt(3.);
 594     for(Int_t j=1;j<=h->GetNbinsX();j++){
 595       h->SetBinContent(j,fhEnvelopeMaxRatio->GetBinContent(j)/sqrtThree-1./sqrtThree);
 596     }
 597     if(fSystUseRMSfromFlatDistr==1)return h;
 598
 599     if(fSystUseRMSfromFlatDistr==2){
 600       TH1D *hRefl=ReflectHisto(h);
 601       TH1D *hBackTo2Pi=DuplicateHistoTo2piRange(hRefl); // In this way we symmetrized to 0-pi
 602       delete h;
 603       delete hRefl;
 604       hBackTo2Pi->SetName("hRelSystFeedDownMin");
 605       return hBackTo2Pi;
 606     }
 607
 608     // case fSystUseRMSfromFlatDistr==3 || ==4
 609     // Do smoothing and then symmetrize
 610     // For the smoothing, assume linear trend every 3 points, starting from transverse region
 611     Int_t jstart=h->FindBin(TMath::Pi()*0.5)-1;
 612     TF1 *f;
 613     if(fSystUseRMSfromFlatDistr==3){
 614       f=new TF1("mypol0","[0]",-TMath::Pi()*0.5,TMath::Pi()*1.5);
 615     }
 616     else if(fSystUseRMSfromFlatDistr==4){
 617       f=new TF1("mypol1","[0]+x*[1]",-TMath::Pi()*0.5,TMath::Pi()*1.5);
 618     }
 619     TH1D *hOut=(TH1D*)h->Clone("hOut");
 620     for(Int_t j=jstart;j>=1;j--){// going towards 0
 621
 622       if(j>=3){
 623         Double_t xmax=h->GetBinLowEdge(j+1)+h->GetBinWidth(j+1);
 624         Double_t xmin=h->GetBinLowEdge(j-1);
 625         f->SetParameter(0,h->GetBinContent(j));
 626         if(fSystUseRMSfromFlatDistr==4)f->SetParameter(1,(h->GetBinContent(j)-h->GetBinContent(j-1))/(h->GetBinCenter(j+1)-h->GetBinCenter(j-1)));
 627         h->Fit(f,"REM","N",xmin,xmax);
 628         hOut->SetBinContent(j,f->Eval(hOut->GetBinCenter(j)));
 629       }
 630       else if(j==2){
 631         Double_t xmax=h->GetBinLowEdge(3)+h->GetBinWidth(3);
 632         Double_t xmin=h->GetBinLowEdge(1);
 633         Double_t y1=h->GetBinContent(1);
 634         Double_t y2=h->GetBinContent(2);
 635         Double_t y3=h->GetBinContent(3);
 636         h->SetBinContent(3,y2);
 637         h->SetBinContent(2,y1);
 638         h->SetBinContent(1,h->GetBinContent(h->GetNbinsX()));
 639         f->SetParameter(0,h->GetBinContent(2));
 640         if(fSystUseRMSfromFlatDistr==4)f->SetParameter(1,(h->GetBinContent(3)-h->GetBinContent(1))/(h->GetBinCenter(3)-h->GetBinCenter(1)));
 641         h->Fit(f,"RLEM","N",xmin,xmax);
 642         hOut->SetBinContent(j,f->Eval(h->GetBinCenter(2)));
 643         h->SetBinContent(3,y3);
 644         h->SetBinContent(2,y2);
 645         h->SetBinContent(1,y1);
 646       }
 647       else if(j==1){
 648         Double_t xmax=h->GetBinLowEdge(3)+h->GetBinWidth(3);
 649         Double_t xmin=h->GetBinLowEdge(1);
 650         Double_t y1=h->GetBinContent(1);
 651         Double_t y2=h->GetBinContent(2);
 652         Double_t y3=h->GetBinContent(3);
 653         h->SetBinContent(3,y1);
 654         h->SetBinContent(2,h->GetBinContent(h->GetNbinsX()));
 655         h->SetBinContent(1,h->GetBinContent(h->GetNbinsX()-1));
 656         f->SetParameter(0,h->GetBinContent(2));
 657         if(fSystUseRMSfromFlatDistr==4) f->SetParameter(1,(h->GetBinContent(3)-h->GetBinContent(1))/(h->GetBinCenter(3)-h->GetBinCenter(1)));
 658         h->Fit(f,"RLEM","N",xmin,xmax);
 659         hOut->SetBinContent(j,f->Eval(h->GetBinCenter(2)));
 660         h->SetBinContent(3,y3);
 661         h->SetBinContent(2,y2);
 662         h->SetBinContent(1,y1);
 663       }
 664
 665     }
 666     for(Int_t j=jstart+1;j<=hOut->GetNbinsX();j++){// now going towards 2pi
 667
 668       if(j<=hOut->GetNbinsX()-2){
 669         Double_t xmax=h->GetBinLowEdge(j+1)+h->GetBinWidth(j+1);
 670         Double_t xmin=h->GetBinLowEdge(j-1);
 671         f->SetParameter(0,h->GetBinContent(j));
 672         if(fSystUseRMSfromFlatDistr==4)f->SetParameter(1,(h->GetBinContent(j)-h->GetBinContent(j-1))/(h->GetBinCenter(j+1)-h->GetBinCenter(j-1)));
 673         h->Fit(f,"RLEM","N",xmin,xmax);
 674         hOut->SetBinContent(j,f->Eval(hOut->GetBinCenter(j)));
 675       }
 676       else if(j==hOut->GetNbinsX()-1){
 677         Double_t xmax=h->GetBinLowEdge(hOut->GetNbinsX())+h->GetBinWidth(hOut->GetNbinsX());
 678         Double_t xmin=h->GetBinLowEdge(hOut->GetNbinsX()-2);
 679         Double_t y1=h->GetBinContent(hOut->GetNbinsX()-2);
 680         Double_t y2=h->GetBinContent(hOut->GetNbinsX()-1);
 681         Double_t y3=h->GetBinContent(hOut->GetNbinsX());
 682         h->SetBinContent(hOut->GetNbinsX()-2,y2);
 683         h->SetBinContent(hOut->GetNbinsX()-1,y3);
 684         h->SetBinContent(hOut->GetNbinsX(),h->GetBinContent(1));
 685         f->SetParameter(0,h->GetBinContent(hOut->GetNbinsX()-1));
 686         if(fSystUseRMSfromFlatDistr==4)f->SetParameter(1,(h->GetBinContent(hOut->GetNbinsX())-h->GetBinContent(hOut->GetNbinsX()-2))/(h->GetBinCenter(hOut->GetNbinsX())-h->GetBinCenter(hOut->GetNbinsX()-2)));
 687         h->Fit(f,"RLEM","N",xmin,xmax);
 688         hOut->SetBinContent(j,f->Eval(h->GetBinCenter(hOut->GetNbinsX()-1)));
 689         h->SetBinContent(hOut->GetNbinsX()-2,y1);
 690         h->SetBinContent(hOut->GetNbinsX()-1,y2);
 691         h->SetBinContent(hOut->GetNbinsX(),y3);
 692       }
 693       else if(j==hOut->GetNbinsX()){
 694         Double_t xmax=h->GetBinLowEdge(hOut->GetNbinsX())+h->GetBinWidth(hOut->GetNbinsX());
 695         Double_t xmin=h->GetBinLowEdge(hOut->GetNbinsX()-2);
 696         Double_t y1=h->GetBinContent(hOut->GetNbinsX()-2);
 697         Double_t y2=h->GetBinContent(hOut->GetNbinsX()-1);
 698         Double_t y3=h->GetBinContent(hOut->GetNbinsX());
 699         h->SetBinContent(hOut->GetNbinsX()-2,y3);
 700         h->SetBinContent(hOut->GetNbinsX()-1,h->GetBinContent(1));
 701         h->SetBinContent(hOut->GetNbinsX(),h->GetBinContent(2));
 702         f->SetParameter(0,h->GetBinContent(hOut->GetNbinsX()-1));
 703         if(fSystUseRMSfromFlatDistr==4) f->SetParameter(1,(h->GetBinContent(hOut->GetNbinsX())-h->GetBinContent(hOut->GetNbinsX()-2))/(h->GetBinCenter(hOut->GetNbinsX())-h->GetBinCenter(hOut->GetNbinsX()-2)));
 704         h->Fit(f,"RLEM","N",xmin,xmax);
 705         hOut->SetBinContent(j,f->Eval(h->GetBinCenter(hOut->GetNbinsX()-1)));
 706         h->SetBinContent(hOut->GetNbinsX()-2,y1);
 707         h->SetBinContent(hOut->GetNbinsX()-1,y2);
 708         h->SetBinContent(hOut->GetNbinsX(),y3);
 709       }
 710
 711     }
 712     // Smoothing done, now simmetrize
 713     TH1D *hRefl=ReflectHisto(hOut);
 714     TH1D *hBackTo2Pi=DuplicateHistoTo2piRange(hRefl); // In this way we symmetrized to 0-pi
 715     delete h;
 716     delete hRefl;
 717     delete hOut;
 718     delete f;
 719     hBackTo2Pi->SetName("hRelSystFeedDownMin");
 720     return hBackTo2Pi;
 721
 722
 723   }
 724
 725   Printf("AliHFCorrelationFDsubtraction: wrong option for getting syst unc");
 726   return 0x0;
 727
 728
 729
 730 }
 731
 732
 733 void AliHFCorrelationFDsubtraction::CalculateEnvelope(){
 734
 735   if(fcAllRatio)delete fcAllRatio;
 736   if(fhEnvelopeMax)delete fhEnvelopeMax;
 737   if(fhEnvelopeMin)delete fhEnvelopeMin;
 738   if(fgrEnvelope)delete   fgrEnvelope;
 739
 740   for(Int_t j=0;j<fLastTemplAdded;j++){
 741     SubtractFeedDown(fhTemplates[j]);
 742   }
 743
 744   fhEnvelopeMaxRatio=(TH1D*)fhDataCorrected[0]->Clone("fhEnvelopeMaxRatio");
 745   fhEnvelopeMaxRatio->Reset();
 746   for(Int_t j=1;j<=  fhEnvelopeMaxRatio->GetNbinsX();j++)  fhEnvelopeMaxRatio->SetBinContent(j,1);
 747   fhEnvelopeMinRatio=(TH1D*)fhDataCorrected[0]->Clone("fhEnvelopeMinRatio");
 748   fhEnvelopeMinRatio->Reset();
 749   for(Int_t j=1;j<=  fhEnvelopeMinRatio->GetNbinsX();j++)  fhEnvelopeMinRatio->SetBinContent(j,1);
 750   fgrEnvelopeRatio=new TGraphAsymmErrors();
 751
 752   fhEnvelopeMax=(TH1D*)fhDataCorrected[0]->Clone("fhEnvelopeMax");
 753   fhEnvelopeMin=(TH1D*)fhDataCorrected[0]->Clone("fhEnvelopeMin");
 754   fgrEnvelope=new TGraphAsymmErrors();
 755
 756
 757   Color_t col[10]={kGray,kRed,kBlue,kGreen,kOrange,kViolet,kCyan,kYellow,kSpring,kMagenta};
 758   fcAllRatio=new TCanvas("cRatioAll","cRatioAll",700,700);
 759   fcAllRatio->cd();
 760   // templ 0 is the reference
 761   fhRatioAllToCentralTempl[0]=(TH1D*)fhDataCorrected[1]->Clone(Form("%sRatioToCentralPred",fhDataCorrected[1]->GetName()));
 762   fhRatioAllToCentralTempl[0]->Divide(fhDataCorrected[0]);
 763
 764   fhRatioAllToCentralTempl[1]=(TH1D*)fhDataCorrected[2]->Clone(Form("%sRatioToCentralPred",fhDataCorrected[2]->GetName()));
 765   fhRatioAllToCentralTempl[1]->Divide(fhDataCorrected[0]);
 766
 767
 768   for(Int_t jb=1;jb<=fhRatioAllToCentralTempl[0]->GetNbinsX();jb++){
 769
 770     fhRatioAllToCentralTempl[0]->SetBinError(jb,0.001);
 771
 772     if(fhRatioAllToCentralTempl[0]->GetBinContent(jb)>1.){
 773       fhEnvelopeMaxRatio->SetBinContent(jb,fhRatioAllToCentralTempl[0]->GetBinContent(jb));
 774       fhEnvelopeMax->SetBinContent(jb,fhDataCorrected[1]->GetBinContent(jb));
 775     }
 776     if(fhRatioAllToCentralTempl[0]->GetBinContent(jb)<1.){
 777       fhEnvelopeMinRatio->SetBinContent(jb,fhRatioAllToCentralTempl[0]->GetBinContent(jb));
 778       fhEnvelopeMin->SetBinContent(jb,fhDataCorrected[1]->GetBinContent(jb));
 779     }
 780
 781     fhRatioAllToCentralTempl[1]->SetBinError(jb,0.001);
 782     if(fhRatioAllToCentralTempl[1]->GetBinContent(jb)>fhEnvelopeMaxRatio->GetBinContent(jb)){
 783         fhEnvelopeMaxRatio->SetBinContent(jb,fhRatioAllToCentralTempl[1]->GetBinContent(jb));
 784         fhEnvelopeMax->SetBinContent(jb,fhDataCorrected[2]->GetBinContent(jb));
 785       }
 786     if(fhRatioAllToCentralTempl[1]->GetBinContent(jb)<fhEnvelopeMinRatio->GetBinContent(jb)){
 787         fhEnvelopeMinRatio->SetBinContent(jb,fhRatioAllToCentralTempl[1]->GetBinContent(jb));
 788         fhEnvelopeMin->SetBinContent(jb,fhDataCorrected[2]->GetBinContent(jb));
 789     }
 790   }
 791
 792
 793   fhRatioAllToCentralTempl[0]->Draw();
 794   fhRatioAllToCentralTempl[1]->Draw("same");
 795
 796   for(Int_t j=1;j<fCountTempl;j++){
 797     // templ j, central
 798     fhRatioAllToCentralTempl[j*3-1]=(TH1D*)fhDataCorrected[j*3]->Clone(Form("%sRatioToCentralPred",fhDataCorrected[j*3]->GetName()));
 799     fhRatioAllToCentralTempl[j*3-1]->Divide(fhDataCorrected[0]);
 800     // templ j, min and max
 801     fhRatioAllToCentralTempl[j*3]=(TH1D*)fhDataCorrected[j*3+1]->Clone(Form("%sRatioToCentralPred",fhDataCorrected[j*3+1]->GetName()));
 802     fhRatioAllToCentralTempl[j*3]->Divide(fhDataCorrected[0]);
 803
 804     fhRatioAllToCentralTempl[j*3+1]=(TH1D*)fhDataCorrected[j*3+2]->Clone(Form("%sRatioToCentralPred",fhDataCorrected[j*3+2]->GetName()));
 805     fhRatioAllToCentralTempl[j*3+1]->Divide(fhDataCorrected[0]);
 806
 807     for(Int_t jb=1;jb<=fhRatioAllToCentralTempl[j*3-1]->GetNbinsX();jb++){
 808
 809       fhRatioAllToCentralTempl[j*3-1]->SetBinError(jb,0.001);
 810       if(fhRatioAllToCentralTempl[j*3-1]->GetBinContent(jb)>fhEnvelopeMaxRatio->GetBinContent(jb)){
 811         fhEnvelopeMaxRatio->SetBinContent(jb,fhRatioAllToCentralTempl[j*3-1]->GetBinContent(jb));
 812         fhEnvelopeMax->SetBinContent(jb,fhDataCorrected[j*3-1]->GetBinContent(jb));
 813       }
 814       if(fhRatioAllToCentralTempl[j*3-1]->GetBinContent(jb)<fhEnvelopeMinRatio->GetBinContent(jb)){
 815         fhEnvelopeMinRatio->SetBinContent(jb,fhRatioAllToCentralTempl[j*3-1]->GetBinContent(jb));
 816         fhEnvelopeMin->SetBinContent(jb,fhDataCorrected[j*3-1]->GetBinContent(jb));
 817       }
 818
 819       fhRatioAllToCentralTempl[j*3]->SetBinError(jb,0.001);
 820       if(fhRatioAllToCentralTempl[j*3]->GetBinContent(jb)>fhEnvelopeMaxRatio->GetBinContent(jb)){
 821         fhEnvelopeMaxRatio->SetBinContent(jb,fhRatioAllToCentralTempl[j*3]->GetBinContent(jb));
 822         fhEnvelopeMax->SetBinContent(jb,fhDataCorrected[j*3]->GetBinContent(jb));
 823       }
 824       if(fhRatioAllToCentralTempl[j*3]->GetBinContent(jb)<fhEnvelopeMinRatio->GetBinContent(jb)){
 825         fhEnvelopeMinRatio->SetBinContent(jb,fhRatioAllToCentralTempl[j*3]->GetBinContent(jb));
 826         fhEnvelopeMin->SetBinContent(jb,fhDataCorrected[j*3]->GetBinContent(jb));
 827       }
 828
 829       fhRatioAllToCentralTempl[j*3+1]->SetBinError(jb,0.001);
 830       if(fhRatioAllToCentralTempl[j*3+1]->GetBinContent(jb)>fhEnvelopeMaxRatio->GetBinContent(jb)){
 831         fhEnvelopeMaxRatio->SetBinContent(jb,fhRatioAllToCentralTempl[j*3+1]->GetBinContent(jb));
 832         fhEnvelopeMax->SetBinContent(jb,fhDataCorrected[j*3+1]->GetBinContent(jb));
 833       }
 834       if(fhRatioAllToCentralTempl[j*3+1]->GetBinContent(jb)<fhEnvelopeMinRatio->GetBinContent(jb)){
 835         fhEnvelopeMinRatio->SetBinContent(jb,fhRatioAllToCentralTempl[j*3+1]->GetBinContent(jb));
 836         fhEnvelopeMin->SetBinContent(jb,fhDataCorrected[j*3+1]->GetBinContent(jb));
 837       }
 838
 839
 840     }
 841
 842     if(j<=10){
 843       fhRatioAllToCentralTempl[j*3-1]->SetLineColor(col[j-1]);
 844       fhRatioAllToCentralTempl[j*3]->SetLineColor(col[j-1]+1);
 845       fhRatioAllToCentralTempl[j*3+1]->SetLineColor(col[j-1]+2);
 846     }
 847     else if(j<=20){
 848       fhRatioAllToCentralTempl[j*3-1]->SetLineColor(col[j-10]+3);
 849       fhRatioAllToCentralTempl[j*3]->SetLineColor(col[j-10]+4);
 850       fhRatioAllToCentralTempl[j*3+1]->SetLineColor(col[j-10]+5);
 851     }
 852     else{
 853       Printf("Too many templates, add more colors, aborting");
 854       return;
 855     }
 856     fhRatioAllToCentralTempl[j*3-1]->Draw("same");
 857     fhRatioAllToCentralTempl[j*3]->Draw("same");
 858     fhRatioAllToCentralTempl[j*3+1]->Draw("same");
 859   }
 860
 861   for(Int_t j=1;j<=fhEnvelopeMin->GetNbinsX();j++){
 862     fhEnvelopeMin->SetBinError(j,0);
 863     fhEnvelopeMinRatio->SetBinError(j,0);
 864     fhEnvelopeMax->SetBinError(j,0);
 865     fhEnvelopeMaxRatio->SetBinError(j,0);
 866   }
 867
 868   fgrEnvelope=GetUncGraphFromHistos(fhDataCorrected[0],fhEnvelopeMin,fhEnvelopeMax);
 869   fgrEnvelopeRatio=GetUncGraphFromHistos(0x0,fhEnvelopeMinRatio,fhEnvelopeMaxRatio);
 870
 871
 872 }
 873
 874
 875 // void AliHFCorrelationFDsubtraction::GetTemplateFromFit(TH1D *h,TH1D *hOut,TString strCanv="cFit",Int_t methodFD=0){
 876
 877 //   // Fit Template from B
 878 //   TCanvas *cFit=new TCanvas(strCanv.Data(),strCanv.Data(),700,700);
 879 //   cFit->cd();
 880 //   h->Draw();
 881
 882 //   if(methodFD==0){
 883 //     TF1 *fitFunction=FitPlotsShort(h,2,0,0);
 884 //     for(Int_t j=1;j<=h->GetNbinsX();j++){
 885 //       hOut->SetBinContent(j,fitFunction->Eval(hOut->GetBinCenter(j)));
 886 //       hOut->SetBinError(j,0);
 887 //     }
 888 //     hOut->SetLineColor(kViolet);
 889 //     hOut->Draw("same");
 890
 891 //   }
 892 //   else if(methodFD==1){
 893
 894
 895 //     Double_t nsybc, ensybc,asybc, easybc;
 896 //     Double_t nsybcMC, ensybcMC,asybcMC, easybcMC;
 897 //     cFit->cd();
 898 //     TF1 *fitFunctionMC=FitPlots(h,2,2,3,nsybcMC, ensybcMC,asybcMC, easybcMC);
 899 //     Double_t baseMC=fitFunctionMC->GetParameter(0);
 900
 901 //     TCanvas *cFitData=new TCanvas(Form("Data%s",strCanv.Data()),"cGetTemplatefromFit_FitData",700,700);
 902 //     cFitData->cd();
 903
 904 //     TH1D *hDataForFit=(TH1D*)hOut->Clone("hDataForFit");
 905 //     hDataForFit->Draw("same");
 906
 907 //     TF1 *fitFunctionData=FitPlots(hDataForFit,1,2,3,nsybc, ensybc,asybc, easybc);
 908 //     Double_t baseData=fitFunctionData->GetParameter(0);
 909
 910 //     for(Int_t j=1;j<=h->GetNbinsX();j++){
 911 //       hOut->SetBinContent(j,fitFunctionMC->Eval(hOut->GetBinCenter(j))+baseData-baseMC);
 912 //       hOut->SetBinError(j,0);
 913 //     }
 914
 915 //     cFit->cd();
 916 //     hOut->SetLineColor(kViolet);
 917 //     hOut->DrawClone("same");
 918
 919 //   }
 920
 921 //   return;
 922 // }
 923
 924
 925
 926 // TH1D* AliHFCorrelationFDsubtraction::SubtracFDrough(Double_t ptmin,Double_t ptmax,Int_t methodSubtr=0,Int_t methodFD=0,Int_t rebin=1,Int_t testAssYieldExt=0,TString correlationDataFile="/Users/administrator/ALICE/CHARM/HFCJ/DCorrelations_Test/2013Jul1PFapprovalNoPrelim/FabioD0/1D_Signal_WithEMCorr_Normal_Charg_OriginSuper_Integrated_Bins10to11_Limits_2_4_TreshPt_0.3_Displ_0.00_Data.root", TString correlationFDtemplFile="/Users/administrator/ALICE/CHARM/HFCJ/DCorrelations_Test/MCtemplate/2013June14Fabio/Plots/1D_Signal_WithEMCorr_MCSample_Normal_Charg_OriginSuper_Integrated_Bins10to11_Limits_2_4_TreshPt_0.3_Displ_0.00_Kine.root",TString spectraMacroOutput="/Users/administrator/ALICE/CHARM/HFCJ/DCorrelations_Test/FeedDownSubtraction/ThirdTest/2013Jun17FabioD0properEff/Dati_CorrectEff/HFPtSpectrum_Nb.root"){
 927
 928 //   /*TString correlationDataFile="/Users/administrator/ALICE/CHARM/HFCJ/DCorrelations_Test/FeedDownSubtraction/SecondTest/PlotCorrelDeta1/1D_Signal_WithEMCorr_Normal_Charg_OriginSuper_Integrated_Bins10to11_Limits_2_4_TreshPt_0.3_Displ_0.00_Data.root",TString correlationFDtemplFile="/Users/administrator/ALICE/CHARM/HFCJ/DCorrelations_Test/MCtemplate/2013June17Fabio/Plots/0.3/1D_Signal_WithEMCorr_MCSample_Normal_Charg_OriginSuper_Integrated_Bins10to11_Limits_2_4_TreshPt_0.3_Displ_0.00_Kine.root",TString spectraMacroOutput="../HFPtSpectrum_Nb.root"){
 929 //    */
 930
 931 //   gROOT->LoadMacro("/Users/administrator/ALICE/CHARM/HFCJ/Macro/FitPlots.C");
 932
 933 //   TFile *fDataCorr=TFile::Open(correlationDataFile.Data(),"READ");
 934 //   TCanvas *cData=fDataCorr->Get("c3");
 935 //   TH1D *hData=(TH1D*)cData->FindObject("hsubtract_norm");
 936 //   hData->Sumw2();
 937 //   hData->Rebin(rebin);
 938 //   hData->Scale(1./(Double_t)rebin);
 939
 940 //   TFile *fFDtemplCorr=TFile::Open(correlationFDtemplFile.Data(),"READ");
 941 //   TCanvas *cFDtempl=fFDtemplCorr->Get("c3");
 942 //   TH1D *hFDtemplFile=(TH1D*)cFDtempl->FindObject("hsubtract_norm_2");
 943 //   TH1D *hPromptTempl=(TH1D*)cFDtempl->FindObject("hsubtract_norm_1");
 944 //   TH1D *hFDtempl;
 945 //   if(methodSubtr==0)hFDtempl=(TH1D*)hFDtemplFile->Clone("hTemplDfromB");
 946 //   else if(methodSubtr==1){
 947 //     hFDtempl=(TH1D*)hData->Clone("hTemplDfromBFitFunc");
 948 //     GetTemplateFromFit(hFDtemplFile,hFDtempl,"cFitFD");
 949 //   }
 950 //   else if(TMath::Abs(methodSubtr)==2){// Increase the baseline of MC template in order to match the baseline observed in data
 951 //     hFDtempl=(TH1D*)hData->Clone("hTemplDfromBFitFunc");
 952 //     GetTemplateFromFit(hFDtemplFile,hFDtempl,"cFitFD",1);
 953 //   }
 954
 955 //   cData->Draw();
 956 //   hFDtempl->SetLineColor(kBlue);
 957 //   hFDtempl->SetMarkerColor(kBlue);
 958 //   hFDtempl->Draw("same");
 959
 960
 961 //   TFile *fSpectrum=TFile::Open(spectraMacroOutput.Data(),"READ");
 962 //   TGraphAsymmErrors *gr=(TGraphAsymmErrors*)fSpectrum->Get("gFcConservative");
 963 //   //   TAxis *ax=gr->GetXaxis();
 964 //   //   Int_t binmin=ax->FindBin(ptmin*1.0001);
 965 //   //   Int_t binmax=ax->FindBin(ptmax*0.9999);
 966
 967
 968
 969 //   Double_t *xgr=gr->GetX();
 970 //   Double_t *elxgr=gr->GetEXlow();
 971 //   Double_t *ehxgr=gr->GetEXhigh();
 972 //   Double_t *fprompt=gr->GetY();
 973 //   Double_t *elfprompt=gr->GetEYlow();
 974 //   Double_t *ehfprompt=gr->GetEYhigh();
 975
 976 //   Int_t binmin=-1,binmax=-1;
 977 //   // Look for right bins
 978 //   for(Int_t j=0;j<gr->GetN();j++){
 979 //     if(binmin<0&&ptmin*1.00001>(xgr[j]-elxgr[j])&&ptmin<(xgr[j]+ehxgr[j]))binmin=j;
 980 //     if(binmax<0&&ptmax>(xgr[j]-elxgr[j])&&ptmax*0.9999<(xgr[j]+ehxgr[j]))binmax=j;
 981
 982 //   }
 983
 984 //   if(TMath::Abs(xgr[binmin]-elxgr[binmin]-ptmin)>0.0001){
 985 //     printf("Bin mismatch: xmin=%f \n",xgr[binmin]-elxgr[binmin]);
 986
 987 //     return;
 988 //   }
 989
 990
 991 //   // DEBUG BIAS:
 992 //   //  fprompt[binmin]=1.;
 993 //   //  elfprompt[binmin]=0.;
 994 //   //  ehfprompt[binmin]=0.;
 995
 996 //   //  hFDtempl->Scale(1./5.);
 997
 998 //   TH1D *hFDCentr=(TH1D*)hFDtempl->Clone("hFDtemplCentr");
 999 //   hFDCentr->Scale(1.-fprompt[binmin]);
1000 //   printf("Centr fsec=%f \n",1.-fprompt[binmin]);
1001
1002 //   TH1D *hDataCentr=(TH1D*)hData->Clone("hDataCentr");
1003 //   hDataCentr->Add(hFDCentr,-1.);// Treatment of Stat unc.
1004 //   hDataCentr->Scale(1./fprompt[binmin]);
1005 //   hDataCentr->SetTitle("central f_{prompt}");
1006
1007 //   TH1D *hFDmax=(TH1D*)hFDtempl->Clone("hFDtemplMax");
1008 //   hFDmax->Scale(1.-(fprompt[binmin]-elfprompt[binmin]));
1009 //   printf("Max fsec=%f \n",1.-(fprompt[binmin]-elfprompt[binmin]));
1010 //   TH1D *hDataMaxFD=(TH1D*)hData->Clone("hDataMaxFD");
1011 //   hDataMaxFD->Add(hFDmax,-1.);// Treatment of Stat unc.
1012 //   hDataMaxFD->Scale(1./(fprompt[binmin]-elfprompt[binmin]));
1013 //   hDataMaxFD->SetTitle("max f_{prompt}");
1014
1015 //   TH1D *hFDmin=(TH1D*)hFDtempl->Clone("hFDtemplMin");
1016 //   hFDmin->Scale(1.-(fprompt[binmin]+ehfprompt[binmin]));
1017 //   printf("Min fsec=%f \n",1.-(fprompt[binmin]+ehfprompt[binmin]));
1018 //   TH1D *hDataMinFD=(TH1D*)hData->Clone("hDataMinFD");
1019 //   hDataMinFD->Add(hFDmin,-1.);// Treatment of Stat unc.
1020 //   hDataMinFD->Scale(1./(fprompt[binmin]+ehfprompt[binmin]));
1021 //   hDataMinFD->SetTitle("min f_{prompt}");
1022
1023
1024
1025 //   TCanvas *cCompare=new TCanvas("cCompare","cCompare",700,700);
1026 //   cCompare->cd();
1027 //   hData->Draw();
1028 //   hData->SetTitle("Uncorrected");
1029 //   hData->SetYTitle("#frac{1}{N_{trig}}#frac{dN}{d#Delta#varphi} (rad^{-1})");
1030 //   hData->SetXTitle("#varphi_{ass}-#varphi_{trig} (rad)");
1031 //   hData->SetMarkerStyle(25);
1032 //   hData->SetMarkerColor(kBlack);
1033 //   hData->SetMarkerSize(1);
1034 //   hData->SetLineColor(kBlack);
1035
1036 //   hDataCentr->SetLineColor(kRed);
1037 //   hDataCentr->SetYTitle("#frac{1}{N_{trig}}#frac{dN}{d#Delta#varphi} (rad^{-1})");
1038 //   hDataCentr->SetXTitle("#varphi_{ass}-#varphi_{trig} (rad)");
1039 //   hDataCentr->SetLineWidth(2);
1040 //   hDataCentr->SetMarkerColor(kRed);
1041 //   hDataCentr->SetMarkerStyle(21);
1042 //   hDataCentr->SetMarkerSize(1);
1043 //   hDataCentr->Draw("sames");
1044
1045 //   hDataMinFD->SetLineColor(kBlue);
1046 //   hDataMinFD->SetMarkerColor(kBlue);
1047 //   hDataMinFD->SetMarkerStyle(23);
1048 //   hDataMinFD->SetMarkerSize(1);
1049 //   hDataMinFD->Draw("sames");
1050
1051 //   hDataMaxFD->SetLineColor(kGreen-6);
1052 //   hDataMaxFD->SetMarkerColor(kGreen-6);
1053 //   hDataMaxFD->SetMarkerStyle(22);
1054 //   hDataMaxFD->SetMarkerSize(1);
1055 //   hDataMaxFD->Draw("sames");
1056
1057 //   TCanvas *cRatio=new TCanvas("cRatio","cRatio",700,700);
1058 //   cRatio->cd();
1059 //   TH1D *hRatioMinToCentr=(TH1D*)hDataMinFD->Clone("hRatioMinToCentr");
1060 //   hRatioMinToCentr->Divide(hDataCentr);
1061 //   for(Int_t j=1;j<=hRatioMinToCentr->GetNbinsX();j++){
1062 //     hRatioMinToCentr->SetBinError(j,0.001);
1063 //   }
1064 //   TH1D *hRatioMaxToCentr=(TH1D*)hDataMaxFD->Clone("hRatioMaxToCentr");
1065 //   hRatioMaxToCentr->Divide(hDataCentr);
1066 //   for(Int_t j=1;j<=hRatioMaxToCentr->GetNbinsX();j++){
1067 //     hRatioMaxToCentr->SetBinError(j,0.001);
1068 //   }
1069
1070 //   hRatioMinToCentr->Draw();
1071 //   hRatioMaxToCentr->Draw("same");
1072
1073
1074
1075 //   TCanvas *cCompareToMC=new TCanvas("cCompareToMC","cCompareToMC",700,700);
1076 //   cCompareToMC->cd();
1077 //   TH1D *hDataCentrCompare=hDataCentr->DrawClone();
1078 //   hDataCentrCompare->SetLineColor(kBlack);
1079 //   hDataCentrCompare->SetMarkerColor(kBlack);
1080 //   hPromptTempl->Draw("same");
1081
1082 //   TCanvas *cCompareToMCFit=new TCanvas("cCompareToMCFit","cCompareToMCFit",700,700);
1083 //   cCompareToMCFit->cd();
1084 //   TH1D *hDataCentrCompare2=(TH1D*)hDataCentr->DrawClone();
1085 //   hDataCentrCompare2->SetLineColor(kBlack);
1086 //   hDataCentrCompare2->SetMarkerColor(kBlack);
1087
1088 //   TH1D *hPromptTemplFit=(TH1D*)hDataCentr->Clone("hPromptDtemplatFit");
1089 //   GetTemplateFromFit(hPromptTempl,hPromptTemplFit,"cFitPrompt",1);
1090 //   cCompareToMCFit->cd();
1091 //   hPromptTemplFit->Draw("same");
1092
1093 //   TCanvas *cRatioToMCFit=new TCanvas("cRatioToMCFit","cRatioToMCFit",700,700);
1094 //   cRatioToMCFit->cd();
1095 //   TH1D *hDataCentrRatio=(TH1D*)hDataCentr->Clone("hDataCentrRatio");
1096 //   hDataCentrRatio->SetLineColor(kBlack);
1097 //   hDataCentrRatio->SetMarkerColor(kBlack);
1098 //   hDataCentrRatio->Divide(hPromptTemplFit);
1099 //   hDataCentrRatio->Draw();
1100
1101 //   TCanvas *cCompareFitToFit=new TCanvas("cCompareFitToFit","cCompareFitToFit",700,700);
1102 //   cCompareFitToFit->cd();
1103 //   TH1D *hDataCentrCompare3=(TH1D*)hDataCentr->DrawClone();
1104 //   hDataCentrCompare3->SetLineColor(kBlack);
1105 //   hDataCentrCompare3->SetMarkerColor(kBlack);
1106
1107 //   Double_t nsybcSt=0.,ensybcSt=0.,asybcSt=0.,easybcSt=0.;
1108 //   Double_t nsySt=0., ensySt=0.,asySt=0., easySt=0.,sigmaNSSt=0.,esigmaNSSt=0.;
1109 //   TF1 *fitSt=FitPlots(hDataCentrCompare2,1,0,0,   nsybcSt, ensybcSt,asybcSt, easybcSt);
1110 //   nsySt=fitSt->GetParameter(1);
1111 //   ensySt=fitSt->GetParError(1);
1112 //   asySt=fitSt->GetParameter(4);
1113 //   easySt=fitSt->GetParError(4);
1114 //   sigmaNSSt=fitSt->GetParameter(3);
1115 //   esigmaNSSt=fitSt->GetParameter(3);
1116
1117 //   hPromptTempl->Draw("sames");
1118
1119 //   if(!testAssYieldExt)  return hDataCentr;
1120
1121 //   // NOW SECTION ON ASSOCIATED YIELD EXTRACTION
1122 //   Double_t nsybc=0., ensybc=0.,asybc=0.,easybc=0.;
1123 //   Double_t nsy=0., ensy=0.,asy=0., easy=0.,sigmaNS=0.,esigmaNS=0.;
1124
1125 //   TH1D *hResidualNSyield=new TH1D("hResidualNSyield","hResidualNSyield",200,0.,10.);
1126 //   TH1D *hResidualASyield=new TH1D("hResidualASyield","hResidualASyield",200,0.,10.);
1127 //   TH1D *hResidualNSsigma=new TH1D("hResidualNSsigma","hResidualNSsigma",200,0.,4.);
1128 //   TH1D *hResidualASsigma=new TH1D("hResidualASsigma","hResidualASsigma",200,0.,4.);
1129
1130 //   TH1D *hResidualNSyieldBC=new TH1D("hResidualNSyieldBC","hResidualNSyieldBC",200,0.,10.);
1131 //   TH1D *hResidualASyieldBC=new TH1D("hResidualASyieldBC","hResidualASyieldBC",200,0.,10.);
1132
1133 //   TCanvas *cResiduals=new TCanvas("cResiduals","cResiduals",700,700);
1134 //   cResiduals->Divide(2,2);
1135 //   cResiduals->cd(1);
1136
1137 //   hResidualNSyield->Draw();
1138 //   hResidualNSyieldBC->SetLineColor(kRed);
1139 //   hResidualNSyieldBC->Draw("same");
1140
1141 //   cResiduals->cd(2);
1142
1143 //   hResidualASyield->Draw();
1144 //   hResidualASyieldBC->SetLineColor(kRed);
1145 //   hResidualASyieldBC->Draw("same");
1146
1147 //   cResiduals->cd(3);
1148 //   hResidualNSsigma->Draw();
1149
1150 //   // START TESTS
1151
1152 //   // FIRST VALUES: STANDARD
1153 //   hResidualNSyield->Fill(nsybcSt);
1154 //   hResidualASyield->Fill(asybcSt);
1155 //   hResidualNSyield->Fill(nsySt);
1156 //   hResidualASyield->Fill(asySt);
1157 //   hResidualNSsigma->Fill(sigmaNSSt);
1158
1159
1160 //   hResidualNSyieldBC->Fill(nsybcSt);
1161 //   hResidualASyieldBC->Fill(asybcSt);
1162
1163
1164 //   // FIX  Baseline at histo min
1165 //   TCanvas *cYieldExtrFixBaselineMin=new TCanvas("cYieldExtrFixBaselineMin","cYieldExtrFixBaselineMin",700,700);
1166 //   cYieldExtrFixBaselineMin->cd();
1167
1168 //   TH1D *hDataCentrCompareBasMin=(TH1D*)hDataCentr->DrawClone();
1169 //   hDataCentrCompareBasMin->SetLineColor(kBlack);
1170 //   hDataCentrCompareBasMin->SetMarkerColor(kBlack);
1171
1172 //   TF1 *fitTest=FitPlots(hDataCentrCompareBasMin,1,1,0,nsybc,ensybc,asybc,easybc);
1173 //   hPromptTempl->Draw("sames");
1174
1175 //   nsy=fitTest->GetParameter(1);
1176 //   ensy=fitTest->GetParError(1);
1177 //   asy=fitTest->GetParameter(4);
1178 //   easy=fitTest->GetParError(4);
1179 //   sigmaNS=fitTest->GetParameter(3);
1180 //   esigmaNS=fitTest->GetParameter(3);
1181
1182 //   hResidualNSyield->Fill(nsybc);
1183 //   hResidualASyield->Fill(asybc);
1184 //   hResidualNSyield->Fill(nsy);
1185 //   hResidualASyield->Fill(asy);
1186 //   hResidualNSsigma->Fill(sigmaNS);
1187
1188 //   hResidualNSyieldBC->Fill(nsybc);
1189 //   hResidualASyieldBC->Fill(asybc);
1190
1191 //   printf("Fix baseline at min ok \n");
1192
1193
1194 //   // FIX  Baseline at average of lower 3 histo points
1195 //   TCanvas *cYieldExtrFixBaseAv3=new TCanvas("cYieldExtrFixBaseAv3","cYieldExtrFixBaseAv3",700,700);
1196 //   cYieldExtrFixBaseAv3->cd();
1197 //   TH1D *hDataCentrCompareBasAv3=(TH1D*)hDataCentr->DrawClone();
1198 //   hDataCentrCompareBasAv3->SetLineColor(kBlack);
1199 //   hDataCentrCompareBasAv3->SetMarkerColor(kBlack);
1200
1201 //   TF1 *fitTest=FitPlots(hDataCentrCompareBasAv3,1,-3,0,nsybc,ensybc,asybc,easybc);
1202 //   hPromptTempl->Draw("sames");
1203
1204 //   nsy=fitTest->GetParameter(1);
1205 //   ensy=fitTest->GetParError(1);
1206 //   asy=fitTest->GetParameter(4);
1207 //   easy=fitTest->GetParError(4);
1208 //   sigmaNS=fitTest->GetParameter(3);
1209 //   esigmaNS=fitTest->GetParameter(3);
1210
1211 //   hResidualNSyield->Fill(nsybc);
1212 //   hResidualASyield->Fill(asybc);
1213 //   hResidualNSyield->Fill(nsy);
1214 //   hResidualASyield->Fill(asy);
1215 //   hResidualNSsigma->Fill(sigmaNS);
1216
1217 //   hResidualNSyieldBC->Fill(nsybc);
1218 //   hResidualASyieldBC->Fill(asybc);
1219
1220 //   printf("Fix baseline at av lower3 ok \n");
1221
1222 //   // FIX  Baseline at average of lower 2 histo points
1223 //   TCanvas *cYieldExtrFixBaseAv2=new TCanvas("cYieldExtrFixBaseAv2","cYieldExtrFixBaseAv2",700,700);
1224 //   cYieldExtrFixBaseAv2->cd();
1225 //   TH1D *hDataCentrCompareBasAv2=(TH1D*)hDataCentr->DrawClone();
1226 //   hDataCentrCompareBasAv2->SetLineColor(kBlack);
1227 //   hDataCentrCompareBasAv2->SetMarkerColor(kBlack);
1228
1229 //   TF1 *fitTest=FitPlots(hDataCentrCompareBasAv2,1,-2,0,nsybc,ensybc,asybc,easybc);
1230 //   hPromptTempl->Draw("sames");
1231
1232 //   nsy=fitTest->GetParameter(1);
1233 //   ensy=fitTest->GetParError(1);
1234 //   asy=fitTest->GetParameter(4);
1235 //   easy=fitTest->GetParError(4);
1236 //   sigmaNS=fitTest->GetParameter(3);
1237 //   esigmaNS=fitTest->GetParameter(3);
1238
1239 //   hResidualNSyield->Fill(nsybc);
1240 //   hResidualASyield->Fill(asybc);
1241 //   hResidualNSyield->Fill(nsy);
1242 //   hResidualASyield->Fill(asy);
1243 //   hResidualNSsigma->Fill(sigmaNS);
1244
1245 //   hResidualNSyieldBC->Fill(nsybc);
1246 //   hResidualASyieldBC->Fill(asybc);
1247
1248
1249 //   printf("Fix baseline at av lower2 ok \n");
1250
1251 //   // FIX  BOTH MEANS
1252 //   TCanvas *cYieldExtrFixBothMeans=new TCanvas("cYieldExtrFixBothMeans","cYieldExtrFixBothMeans",700,700);
1253 //   cYieldExtrFixBothMeans->cd();
1254 //   TH1D *hDataCentrCompareFMB=(TH1D*)hDataCentr->DrawClone();
1255 //   hDataCentrCompareFMB->SetLineColor(kBlack);
1256 //   hDataCentrCompareFMB->SetMarkerColor(kBlack);
1257
1258 //   TF1 *fitTest=FitPlots(hDataCentrCompareFMB,1,0,3,nsybc,ensybc,asybc,easybc);
1259 //   hPromptTempl->Draw("sames");
1260
1261 //   nsy=fitTest->GetParameter(1);
1262 //   ensy=fitTest->GetParError(1);
1263 //   asy=fitTest->GetParameter(4);
1264 //   easy=fitTest->GetParError(4);
1265 //   sigmaNS=fitTest->GetParameter(3);
1266 //   esigmaNS=fitTest->GetParameter(3);
1267
1268 //   hResidualNSyield->Fill(nsybc);
1269 //   hResidualASyield->Fill(asybc);
1270 //   hResidualNSyield->Fill(nsy);
1271 //   hResidualASyield->Fill(asy);
1272 //   hResidualNSsigma->Fill(sigmaNS);
1273
1274 //   hResidualNSyieldBC->Fill(nsybc);
1275 //   hResidualASyieldBC->Fill(asybc);
1276
1277 //   printf("Fix both means ok \n");
1278 //   return;
1279
1280 //   // FIX  NS MEAN
1281 //   TCanvas *cYieldExtrFixNSMean=new TCanvas("cYieldExtrFixNSMean","cYieldExtrFixNSMean",700,700);
1282 //   cYieldExtrFixBothMeans->cd();
1283 //   TH1D *hDataCentrCompareFMNS=(TH1D*)hDataCentr->DrawClone();
1284 //   hDataCentrCompareFMNS->SetLineColor(kBlack);
1285 //   hDataCentrCompareFMNS->SetMarkerColor(kBlack);
1286
1287 //   TF1 *fitTest=FitPlots(hDataCentrCompareFMNS,1,0,1,nsybc,ensybc,asybc,easybc);
1288 //   hPromptTempl->Draw("sames");
1289
1290 //   nsy=fitTest->GetParameter(1);
1291 //   ensy=fitTest->GetParError(1);
1292 //   asy=fitTest->GetParameter(4);
1293 //   easy=fitTest->GetParError(4);
1294 //   sigmaNS=fitTest->GetParameter(3);
1295 //   esigmaNS=fitTest->GetParameter(3);
1296
1297 //   hResidualNSyield->Fill(nsybc);
1298 //   hResidualASyield->Fill(asybc);
1299 //   hResidualNSyield->Fill(nsy);
1300 //   hResidualASyield->Fill(asy);
1301 //   hResidualNSsigma->Fill(sigmaNS);
1302
1303 //   hResidualNSyieldBC->Fill(nsybc);
1304 //   hResidualASyieldBC->Fill(asybc);
1305
1306 //   printf("Fix NS mean ok \n");
1307
1308 //   // FIX  AS MEAN
1309 //   TCanvas *cYieldExtrFixASMean=new TCanvas("cYieldExtrFixASMean","cYieldExtrFixASMean",700,700);
1310 //   cYieldExtrFixBothMeans->cd();
1311 //   TH1D *hDataCentrCompareFMAS=(TH1D*)hDataCentr->DrawClone();
1312 //   hDataCentrCompareFMAS->SetLineColor(kBlack);
1313 //   hDataCentrCompareFMAS->SetMarkerColor(kBlack);
1314
1315 //   TF1 *fitTest=FitPlots(hDataCentrCompareFMAS,1,0,2,nsybc,ensybc,asybc,easybc);
1316 //   hPromptTempl->Draw("sames");
1317
1318 //   nsy=fitTest->GetParameter(1);
1319 //   ensy=fitTest->GetParError(1);
1320 //   asy=fitTest->GetParameter(4);
1321 //   easy=fitTest->GetParError(4);
1322 //   sigmaNS=fitTest->GetParameter(3);
1323 //   esigmaNS=fitTest->GetParameter(3);
1324
1325 //   hResidualNSyield->Fill(nsybc);
1326 //   hResidualASyield->Fill(asybc);
1327 //   hResidualNSyield->Fill(nsy);
1328 //   hResidualASyield->Fill(asy);
1329 //   hResidualNSsigma->Fill(sigmaNS);
1330
1331 //   hResidualNSyieldBC->Fill(nsybc);
1332 //   hResidualASyieldBC->Fill(asybc);
1333
1334 //   printf("Fix AS mean ok \n");
1335
1336
1337 // }