PWGHF/correlationHF/AliHFDmesonCorrAverage.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
  17 /* $Id: $ */
  18
  19 /////////////////////////////////////////////////////////////
  20 // class to average D meson -hadron correlations
  21 //
  22 // Author: A. Rossi, andrea.rossi@cern.ch
  23 /////////////////////////////////////////////////////////////
  24
  25 #include "AliHFDmesonCorrAverage.h"
  26 #include "AliHFDhadronCorrSystUnc.h"
  27 #include <TMath.h>
  28 #include <TH1D.h>
  29 #include <TCanvas.h>
  30 #include <TGraphAsymmErrors.h>
  31
  32
  33 //using std::cout;
  34 //using std::endl;
  35
  36 ClassImp(AliHFDmesonCorrAverage)
  37
  38 AliHFDmesonCorrAverage::AliHFDmesonCorrAverage() :
  39 TNamed(),
  40   fSystDzero(0),
  41   fSystDstar(0),
  42   fSystDplus(0),
  43   fSystDaverage(0),
  44   fincludeDzero(kFALSE),
  45   fincludeDstar(kFALSE),
  46   fincludeDplus(kFALSE),
  47   fmethod(),
  48   fptminD(),
  49   fptmaxD(),
  50   fptminAsso(),
  51   fptmaxAsso(),
  52   fhDzero(0x0),
  53   fhDstar(0x0),
  54   fhDplus(0x0),
  55   fhDaverage(0x0),
  56   fgrTotSystAverage(0x0),
  57   fgrFDSystAverage(0x0),
  58   fgrNonFDSystAverage(0x0),
  59   fweightsDzeroStat(0x0),
  60   fweightsDstarStat(0x0),
  61   fweightsDplusStat(0x0),
  62   fweightsDzeroSystYield(0x0),
  63   fweightsDstarSystYield(0x0),
  64   fweightsDplusSystYield(0x0),
  65   fweightsDzeroSystBkg(0x0),
  66   fweightsDstarSystBkg(0x0),
  67   fweightsDplusSystBkg(0x0),
  68   fnbinsphi(),
  69   fsys(-2),
  70   fyear(-2),
  71   fSystAlreadySet(kFALSE),
  72   fArithmeticAverage(kFALSE),
  73   fhUsedWeightsDzero(0x0),
  74   fhUsedWeightsDstar(0x0),
  75   fhUsedWeightsDplus(0x0)
  76 {// standard constructor
  77
  78 }
  79
  80
  81 AliHFDmesonCorrAverage::AliHFDmesonCorrAverage(const char* name) :
  82   TNamed(name,name),
  83   fSystDzero(0),
  84   fSystDstar(0),
  85   fSystDplus(0),
  86   fSystDaverage(0),
  87   fincludeDzero(kFALSE),
  88   fincludeDstar(kFALSE),
  89   fincludeDplus(kFALSE),
  90   fmethod(),
  91   fptminD(),
  92   fptmaxD(),
  93   fptminAsso(),
  94   fptmaxAsso(),
  95   fhDzero(0x0),
  96   fhDstar(0x0),
  97   fhDplus(0x0),
  98   fhDaverage(0x0),
  99   fgrTotSystAverage(0x0),
 100   fgrFDSystAverage(0x0),
 101   fgrNonFDSystAverage(0x0),
 102   fweightsDzeroStat(0x0),
 103   fweightsDstarStat(0x0),
 104   fweightsDplusStat(0x0),
 105   fweightsDzeroSystYield(0x0),
 106   fweightsDstarSystYield(0x0),
 107   fweightsDplusSystYield(0x0),
 108   fweightsDzeroSystBkg(0x0),
 109   fweightsDstarSystBkg(0x0),
 110   fweightsDplusSystBkg(0x0),
 111   fnbinsphi(),
 112   fsys(-2),
 113   fyear(-2),
 114   fSystAlreadySet(kFALSE),
 115   fArithmeticAverage(kFALSE),
 116   fhUsedWeightsDzero(0x0),
 117   fhUsedWeightsDstar(0x0),
 118   fhUsedWeightsDplus(0x0)
 119 {// default constructor
 120
 121 }
 122
 123 AliHFDmesonCorrAverage::~AliHFDmesonCorrAverage(){
 124
 125   delete fSystDzero;
 126   delete fSystDstar;
 127   delete fSystDplus;
 128   delete fSystDaverage;
 129   delete fhDzero;
 130   delete fhDstar;
 131   delete fhDplus;
 132   delete fhDaverage;
 133   delete fgrTotSystAverage;
 134   delete fgrFDSystAverage;
 135   delete fgrNonFDSystAverage;
 136   delete fweightsDzeroStat;
 137   delete fweightsDstarStat;
 138   delete fweightsDplusStat;
 139   delete fweightsDzeroSystYield;
 140   delete fweightsDstarSystYield;
 141   delete fweightsDplusSystYield;
 142   delete fweightsDzeroSystBkg;
 143   delete fweightsDstarSystBkg;
 144   delete fweightsDplusSystBkg;
 145   delete  fhUsedWeightsDzero;
 146   delete fhUsedWeightsDstar;
 147   delete fhUsedWeightsDplus;
 148
 149 }
 150
 151
 152
 153 Bool_t AliHFDmesonCorrAverage::InitSystematicUncertainty(Int_t system,Int_t year){
 154   if(!fSystAlreadySet){
 155     if(system!=-1){
 156       if(system==0){ //pp
 157         if(fincludeDzero){
 158           if(year==2010){
 159             fSystDzero=new AliHFDhadronCorrSystUnc("fSystDzero");
 160             fSystDzero->InitStandardUncertaintiesPP2010(0,(fptminD+fptmaxD)*0.5,fptminAsso);
 161             fSystDzero->BuildTotalUncHisto();
 162             fSystDzero->BuildTotalNonFDUncHisto();
 163             fSystDzero->BuildGraphsUnc(fhDzero);
 164             fSystDzero->BuildGraphsRelUnc();
 165           }
 166           else {
 167             Printf("No values for syst unc foreseen for this dataset");
 168           }
 169         }
 170
 171         if(fincludeDstar){
 172           if(year==2010){
 173             fSystDstar=new AliHFDhadronCorrSystUnc("fSystDstar");
 174             fSystDstar->InitStandardUncertaintiesPP2010(1,(fptminD+fptmaxD)*0.5,fptminAsso);
 175             fSystDstar->BuildTotalUncHisto();
 176             fSystDstar->BuildTotalNonFDUncHisto();
 177             fSystDstar->BuildGraphsUnc(fhDstar);
 178             fSystDstar->BuildGraphsRelUnc();
 179           }
 180           else {
 181             Printf("No values for syst unc foreseen for this dataset");
 182           }
 183         }
 184
 185         if(fincludeDplus){
 186           if(year==2010){
 187             fSystDplus=new AliHFDhadronCorrSystUnc("fSystDplus");
 188             fSystDplus->InitStandardUncertaintiesPP2010(2,(fptminD+fptmaxD)*0.5,fptminAsso);
 189             fSystDplus->BuildTotalUncHisto();
 190             fSystDplus->BuildTotalNonFDUncHisto();
 191             fSystDplus->BuildGraphsUnc(fhDplus);
 192             fSystDplus->BuildGraphsRelUnc();
 193           }
 194           else {
 195             Printf("No values for syst unc foreseen for this dataset");
 196           }
 197         }
 198       }
 199       else if(system==1){ //p-Pb
 200         if(fincludeDzero){
 201           if(year==2013){
 202             fSystDzero=new AliHFDhadronCorrSystUnc("fSystDzero");
 203             fSystDzero->InitStandardUncertaintiesPPb2013(0,(fptminD+fptmaxD)*0.5,fptminAsso);
 204             fSystDzero->BuildTotalUncHisto();
 205             fSystDzero->BuildTotalNonFDUncHisto();
 206             fSystDzero->BuildGraphsUnc(fhDzero);
 207             fSystDzero->BuildGraphsRelUnc();
 208           }
 209           else {
 210             Printf("No values for syst unc foreseen for this dataset");
 211           }
 212         }
 213         if(fincludeDstar){
 214           if(year==2013){
 215             fSystDstar=new AliHFDhadronCorrSystUnc("fSystDstar");
 216             fSystDstar->InitStandardUncertaintiesPPb2013(1,(fptminD+fptmaxD)*0.5,fptminAsso);
 217             fSystDstar->BuildTotalUncHisto();
 218             fSystDstar->BuildTotalNonFDUncHisto();
 219             fSystDstar->BuildGraphsUnc(fhDstar);
 220             fSystDstar->BuildGraphsRelUnc();
 221           }
 222           else {
 223             Printf("No values for syst unc foreseen for this dataset");
 224           }
 225         }
 226
 227         if(fincludeDplus){
 228           if(year==2013){
 229             fSystDplus=new AliHFDhadronCorrSystUnc("fSystDplus");
 230             fSystDplus->InitStandardUncertaintiesPPb2013(2,(fptminD+fptmaxD)*0.5,fptminAsso);
 231             fSystDplus->BuildTotalUncHisto();
 232             fSystDplus->BuildTotalNonFDUncHisto();
 233             fSystDplus->BuildGraphsUnc(fhDplus);
 234             fSystDplus->BuildGraphsRelUnc();
 235           }
 236           else {
 237             Printf("No values for syst unc foreseen for this dataset");
 238           }
 239         }
 240       }
 241       else {
 242         Printf("Cannot initiate syst uncertainties: wrong system selected");
 243         return kFALSE;
 244       }
 245     }
 246   }
 247
 248   fSystDaverage=new AliHFDhadronCorrSystUnc("fSystDaverage");
 249   if(fincludeDzero)fSystDaverage->SetHistoTemplate(fSystDzero->GetHistoTemplate(),"fhDeltaPhiTemplateDaverage");
 250   else if(fincludeDstar)fSystDaverage->SetHistoTemplate(fSystDstar->GetHistoTemplate(),"fhDeltaPhiTemplateDaverage");
 251   else if(fincludeDplus)fSystDaverage->SetHistoTemplate(fSystDplus->GetHistoTemplate(),"fhDeltaPhiTemplateDaverage");
 252
 253   fSystDaverage->InitEmptyHistosFromTemplate();
 254
 255   return kTRUE;
 256 }
 257
 258
 259 void AliHFDmesonCorrAverage::InitAverageHisto(TH1D *h){
 260   fhDaverage=(TH1D*)h->Clone("fhDaverage");
 261   fhDaverage->SetXTitle("#Delta#varphi = #varphi_{assoc} - #varphi_{D}");
 262   fhDaverage->SetYTitle("#frac{1}{N_{D}}#frac{dN^{h}}{d#Delta#varphi} (rad^{-1})");
 263   fhDaverage->SetTitle("");
 264   fhDaverage->SetMinimum(0);
 265
 266   fhDaverage->Reset(0);
 267   fhDaverage->SetLineWidth(2);
 268   fhDaverage->SetLineColor(kBlack);
 269   fhDaverage->SetMarkerStyle(20);
 270   fhDaverage->SetMarkerSize(1.2);
 271   fhDaverage->SetMarkerColor(kBlack);
 272
 273   // The following histos are created here to use the final binning
 274   fhUsedWeightsDzero=(TH1D*)h->Clone("fhUsedWeightsDzero");
 275   fhUsedWeightsDzero->SetTitle("Dzero final weights used");
 276   fhUsedWeightsDzero->SetXTitle("#Delta#varphi = #varphi_{assoc} - #varphi_{D}");
 277   fhUsedWeightsDzero->SetYTitle("weight");
 278
 279   fhUsedWeightsDplus=(TH1D*)h->Clone("fhUsedWeightsDplus");
 280   fhUsedWeightsDplus->SetTitle("Dplus final weights used");
 281   fhUsedWeightsDplus->SetXTitle("#Delta#varphi = #varphi_{assoc} - #varphi_{D}");
 282   fhUsedWeightsDplus->SetYTitle("weight");
 283
 284   fhUsedWeightsDstar=(TH1D*)h->Clone("fhUsedWeightsDstar");
 285   fhUsedWeightsDstar->SetTitle("Dstar final weights used");
 286   fhUsedWeightsDstar->SetXTitle("#Delta#varphi = #varphi_{assoc} - #varphi_{D}");
 287   fhUsedWeightsDstar->SetYTitle("weight");
 288
 289
 290
 291 }
 292
 293 void AliHFDmesonCorrAverage::CalculateAverage(){
 294   if(fmethod!=10&&fmethod!=20){
 295     Printf("This option for the calculation of the average has not been implemented yet");
 296     return;
 297   }
 298   // check that histos were set
 299   if(fincludeDzero&&(!fhDzero)){
 300     Printf("Dzero histo not set");
 301     return;
 302   }
 303   if(fincludeDstar&&(!fhDstar)){
 304     Printf("Dstar histo not set");
 305     return;
 306   }
 307   if(fincludeDplus&&(!fhDplus)){
 308     Printf("Dplus histo not set");
 309     return;
 310   }
 311
 312
 313   // check that they have the same binning
 314   if(fincludeDzero&&fincludeDstar){
 315     if(fhDzero->GetNbinsX()!=fhDstar->GetNbinsX()){
 316       Printf("Dzero and Dstar histos have different binning");
 317       return ;
 318     }
 319     fnbinsphi=fhDzero->GetNbinsX();
 320     for(Int_t j=1;j<=fnbinsphi;j++){
 321       if(TMath::Abs(fhDzero->GetBinLowEdge(j)-fhDstar->GetBinLowEdge(j))>0.001){
 322         Printf("Dzero and Dstar histos have different binning");
 323         return;
 324       }
 325     }
 326     InitAverageHisto(fhDzero);
 327   }
 328   if(fincludeDzero&&fincludeDplus){
 329     if(fhDzero->GetNbinsX()!=fhDplus->GetNbinsX()){
 330     Printf("Dzero and Dplus histos have different binning");
 331     return ;
 332     }
 333     fnbinsphi=fhDzero->GetNbinsX();
 334     for(Int_t j=1;j<=fnbinsphi;j++){
 335       if(TMath::Abs(fhDzero->GetBinLowEdge(j)-fhDplus->GetBinLowEdge(j))>0.001){
 336         Printf("Dzero and Dplus histos have different binning");
 337         return;
 338       }
 339     }
 340     if(!fhDaverage)InitAverageHisto(fhDzero);
 341   }
 342   if(fincludeDstar&&fincludeDplus){
 343     if(fhDstar->GetNbinsX()!=fhDplus->GetNbinsX()){
 344       Printf("Dstar and Dplus histos have different binning");
 345       return ;
 346     }
 347     fnbinsphi=fhDstar->GetNbinsX();
 348     for(Int_t j=1;j<=fnbinsphi;j++){
 349       if(TMath::Abs(fhDstar->GetBinLowEdge(j)-fhDplus->GetBinLowEdge(j))>0.001){
 350         Printf("Dplus and Dstar histos have different binning");
 351         return;
 352       }
 353     }
 354     if(!fhDaverage)InitAverageHisto(fhDstar);
 355   }
 356
 357   if(!(fincludeDstar||fincludeDplus||fincludeDzero)){
 358     Printf("No mesons choosen for average");
 359     return ;
 360   }
 361   Int_t nmeson=0;
 362   if(fincludeDzero)nmeson++;
 363   if(fincludeDstar)nmeson++;
 364   if(fincludeDplus)nmeson++;
 365
 366   printf("Settin syst \n");
 367
 368   if(!InitSystematicUncertainty(fsys,fyear)){
 369     Printf("Cannot init syst uncertainties \n");
 370     return;
 371   }
 372
 373   SetWeights();
 374   printf("Weights set \n");
 375
 376   for(Int_t j=1;j<=fnbinsphi;j++){
 377     Double_t value=0.,errStatValue=0.,errSystYieldValue=0.,errSystBkgValue=0.,weight=0.;
 378     Double_t systMCclosureMin=0.,systMCcorrectionsMin=0.,systFDmin=0.,systMCDefficiencyMin=0.,systSecContaminationMin=0.;
 379     Double_t systMCclosureMax=0.,systMCcorrectionsMax=0.,systFDmax=0.,systMCDefficiencyMax=0.,systSecContaminationMax=0.;
 380     Double_t sumweights=0;
 381     Double_t errsyst;
 382
 383     if(fincludeDzero){
 384       // stat error + yield unc + bkg subtr
 385       if(fArithmeticAverage){
 386         weight=1./nmeson;
 387       }
 388       else {
 389         weight=1./(1./fweightsDzeroStat[j-1]+1./fweightsDzeroSystYield[j-1]+1./fweightsDzeroSystBkg[j-1]);// need to do this way since we stored separately the stat and syst weight (=1/variance)
 390       }
 391       fhUsedWeightsDzero->SetBinContent(j,weight);
 392
 393       value+=fhDzero->GetBinContent(j)*weight;
 394       errStatValue+=1./fweightsDzeroStat[j-1]*weight*weight;
 395       errSystYieldValue+=1./fweightsDzeroSystYield[j-1]*weight*weight;
 396       errSystBkgValue+=1./fweightsDzeroSystBkg[j-1]*weight*weight;
 397       sumweights+=weight;
 398
 399       printf("Dzero the value is: %f, weight: %f \n",value, weight);
 400       // MCcorrections  : associated tracks, correlated
 401       errsyst=TMath::Abs(fhDzero->GetBinContent(j)*fSystDzero->GetHistoMCcorrectionsMin()->GetBinContent(j));
 402       systMCcorrectionsMin+=errsyst*weight;
 403       errsyst=fhDzero->GetBinContent(j)*fSystDzero->GetHistoMCcorrectionsMax()->GetBinContent(j);
 404       systMCcorrectionsMax+=errsyst*weight;
 405
 406       printf(" Dzero trackeff the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDzero->GetHistoMCcorrectionsMin()->GetBinContent(j), systMCcorrectionsMin);
 407       printf(" Dzero trackeff the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDzero->GetHistoMCcorrectionsMax()->GetBinContent(j), systMCcorrectionsMax);
 408
 409       // MCDefficiency  : uncorrelated
 410       errsyst=fhDzero->GetBinContent(j)*fSystDzero->GetHistoMCDefficiencyMin()->GetBinContent(j);
 411       systMCDefficiencyMin+=errsyst*errsyst*weight*weight;
 412       errsyst=fhDzero->GetBinContent(j)*fSystDzero->GetHistoMCDefficiencyMax()->GetBinContent(j);
 413       systMCDefficiencyMax+=errsyst*errsyst*weight*weight;
 414
 415       // SecContamination  : correlated
 416       errsyst=TMath::Abs(fhDzero->GetBinContent(j)*fSystDzero->GetHistoSecContaminationMin()->GetBinContent(j));
 417       systSecContaminationMin+=errsyst*weight;
 418       errsyst=fhDzero->GetBinContent(j)*fSystDzero->GetHistoSecContaminationMax()->GetBinContent(j);
 419       systSecContaminationMax+=errsyst*weight;
 420
 421       // MC closure: fully correlated
 422       errsyst=TMath::Abs(fhDzero->GetBinContent(j)*fSystDzero->GetHistoMCclosureTestMin()->GetBinContent(j));
 423       systMCclosureMin+=errsyst*weight;
 424       errsyst=fhDzero->GetBinContent(j)*fSystDzero->GetHistoMCclosureTestMax()->GetBinContent(j);
 425       systMCclosureMax+=errsyst*weight;
 426
 427       printf(" Dzero Mcclosure the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDzero->GetHistoMCclosureTestMin()->GetBinContent(j), systMCclosureMin);
 428       printf(" Dzero Mcclosure the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDzero->GetHistoMCclosureTestMax()->GetBinContent(j), systMCclosureMax);
 429
 430       // FD (assume full correlations)
 431       errsyst=TMath::Abs(fhDzero->GetBinContent(j)*fSystDzero->GetHistoFDmin()->GetBinContent(j));
 432       systFDmin+=errsyst*weight;
 433       errsyst=fhDzero->GetBinContent(j)*fSystDzero->GetHistoFDmax()->GetBinContent(j);
 434       systFDmax+=errsyst*weight;
 435       printf(" Dzero feeddown the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDzero->GetHistoFDmin()->GetBinContent(j), systFDmin);
 436       printf(" Dzero feeddown the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDzero->GetHistoFDmax()->GetBinContent(j), systFDmax);
 437
 438       printf("Dzero the value is: %f, weight: %f \n",value, weight);
 439     }
 440     if(fincludeDstar){
 441       if(fArithmeticAverage){
 442         weight=1./nmeson;
 443       }
 444       else{
 445         // stat error + yield unc + bkg subtr
 446         weight=1./(1./fweightsDstarStat[j-1]+1./fweightsDstarSystYield[j-1]+1./fweightsDstarSystBkg[j-1]);// need to do this way since we stored separately the stat and syst weight (=1/variance)
 447       }
 448       fhUsedWeightsDstar->SetBinContent(j,weight);
 449       value+=fhDstar->GetBinContent(j)*weight;
 450       errStatValue+=1./fweightsDstarStat[j-1]*weight*weight;
 451       errSystYieldValue+=1./fweightsDstarSystYield[j-1]*weight*weight;
 452       errSystBkgValue+=1./fweightsDstarSystBkg[j-1]*weight*weight;
 453       sumweights+=weight;
 454
 455       printf("Dstar the value is: %f, weight: %f \n",value, weight);
 456
 457       // MCcorrections: associated tracks, correlated
 458       errsyst=TMath::Abs(fhDstar->GetBinContent(j)*fSystDstar->GetHistoMCcorrectionsMin()->GetBinContent(j));
 459       systMCcorrectionsMin+=errsyst*weight;
 460       errsyst=TMath::Abs(fhDstar->GetBinContent(j)*fSystDstar->GetHistoMCcorrectionsMax()->GetBinContent(j));
 461       systMCcorrectionsMax+=errsyst*weight;
 462
 463       printf(" Dstar trackeff the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDstar->GetHistoMCcorrectionsMin()->GetBinContent(j), systMCcorrectionsMin);
 464       printf(" Dstar trackeff the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDstar->GetHistoMCcorrectionsMax()->GetBinContent(j), systMCcorrectionsMax);
 465
 466       // MCDefficiency  : uncorrelated
 467       errsyst=fhDstar->GetBinContent(j)*fSystDstar->GetHistoMCDefficiencyMin()->GetBinContent(j);
 468       systMCDefficiencyMin+=errsyst*errsyst*weight*weight;
 469       errsyst=fhDstar->GetBinContent(j)*fSystDstar->GetHistoMCDefficiencyMax()->GetBinContent(j);
 470       systMCDefficiencyMax+=errsyst*errsyst*weight*weight;
 471
 472       // SecContamination  : correlated
 473       errsyst=TMath::Abs(fhDstar->GetBinContent(j)*fSystDstar->GetHistoSecContaminationMin()->GetBinContent(j));
 474       systSecContaminationMin+=errsyst*weight;
 475       errsyst=fhDstar->GetBinContent(j)*fSystDstar->GetHistoSecContaminationMax()->GetBinContent(j);
 476       systSecContaminationMax+=errsyst*weight;
 477
 478       // MC closure
 479       errsyst=TMath::Abs(fhDstar->GetBinContent(j)*fSystDstar->GetHistoMCclosureTestMin()->GetBinContent(j));
 480       systMCclosureMin+=errsyst*weight;
 481       errsyst=fhDstar->GetBinContent(j)*fSystDstar->GetHistoMCclosureTestMax()->GetBinContent(j);
 482       systMCclosureMax+=errsyst*weight;
 483
 484       printf(" Dstar Mcclosure the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDstar->GetHistoMCclosureTestMin()->GetBinContent(j), systMCclosureMin);
 485       printf(" Dstar Mcclosure the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDstar->GetHistoMCclosureTestMax()->GetBinContent(j), systMCclosureMax);
 486
 487
 488       // FD (assume full correlations)
 489       errsyst=TMath::Abs(fhDstar->GetBinContent(j)*fSystDstar->GetHistoFDmin()->GetBinContent(j));
 490       systFDmin+=errsyst*weight;
 491       errsyst=fhDstar->GetBinContent(j)*fSystDstar->GetHistoFDmax()->GetBinContent(j);
 492       systFDmax+=errsyst*weight;
 493
 494       printf(" Dstar feeddown the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDstar->GetHistoFDmin()->GetBinContent(j), systFDmin);
 495       printf(" Dstar feeddown the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDstar->GetHistoFDmax()->GetBinContent(j), systFDmax);
 496
 497       printf("Dstar the value is: %f, weight: %f \n",value, weight);
 498     }
 499     if(fincludeDplus){
 500       if(fArithmeticAverage){
 501         weight=1./nmeson;
 502       }
 503       else{
 504         // stat error + yield unc + bkg subtr
 505         weight=1./(1./fweightsDplusStat[j-1]+1./fweightsDplusSystYield[j-1]+1./fweightsDplusSystBkg[j-1]);// need to do this way since we stored separately the stat and syst weight (=1/variance)
 506       }
 507       fhUsedWeightsDplus->SetBinContent(j,weight);
 508       value+=fhDplus->GetBinContent(j)*weight;
 509       errStatValue+=1./fweightsDplusStat[j-1]*weight*weight;
 510       errSystYieldValue+=1./fweightsDplusSystYield[j-1]*weight*weight;
 511       errSystBkgValue+=1./fweightsDplusSystBkg[j-1]*weight*weight;
 512       sumweights+=weight;
 513
 514       printf("Dplus the value is: %f, weight: %f \n",value, weight);
 515
 516       // MCcorrections  : associated tracsk, correlated
 517       errsyst=TMath::Abs(fhDplus->GetBinContent(j)*fSystDplus->GetHistoMCcorrectionsMin()->GetBinContent(j));
 518       systMCcorrectionsMin+=errsyst*weight;
 519       errsyst=fhDplus->GetBinContent(j)*fSystDplus->GetHistoMCcorrectionsMax()->GetBinContent(j);
 520       systMCcorrectionsMax+=errsyst*weight;
 521
 522       printf(" Dplus trackeff the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoMCcorrectionsMin()->GetBinContent(j), systMCcorrectionsMin);
 523       printf(" Dplus trackeff the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoMCcorrectionsMax()->GetBinContent(j), systMCcorrectionsMax);
 524
 525       // MCDefficiency  : uncorrelated
 526       errsyst=fhDplus->GetBinContent(j)*fSystDplus->GetHistoMCDefficiencyMin()->GetBinContent(j);
 527       systMCDefficiencyMin+=errsyst*errsyst*weight*weight;
 528       errsyst=fhDplus->GetBinContent(j)*fSystDplus->GetHistoMCDefficiencyMax()->GetBinContent(j);
 529       systMCDefficiencyMax+=errsyst*errsyst*weight*weight;
 530
 531       printf(" Dplus cutvariat the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoMCDefficiencyMin()->GetBinContent(j), systMCDefficiencyMin);
 532       printf(" Dplus cutvariat the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoMCDefficiencyMax()->GetBinContent(j), systMCDefficiencyMax);
 533
 534       // SecContamination  : correlated
 535       errsyst=TMath::Abs(fhDplus->GetBinContent(j)*fSystDplus->GetHistoSecContaminationMin()->GetBinContent(j));
 536       systSecContaminationMin+=errsyst*weight;
 537       errsyst=fhDplus->GetBinContent(j)*fSystDplus->GetHistoSecContaminationMax()->GetBinContent(j);
 538       systSecContaminationMax+=errsyst*weight;
 539
 540       // MC closure
 541       errsyst=TMath::Abs(fhDplus->GetBinContent(j)*fSystDplus->GetHistoMCclosureTestMin()->GetBinContent(j));
 542       systMCclosureMin+=errsyst*weight;
 543       errsyst=fhDplus->GetBinContent(j)*fSystDplus->GetHistoMCclosureTestMax()->GetBinContent(j);
 544       systMCclosureMax+=errsyst*weight;
 545
 546       printf(" Dplus Mcclosure the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoMCclosureTestMin()->GetBinContent(j), systMCclosureMin);
 547       printf(" Dplus Mcclosure the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoMCclosureTestMax()->GetBinContent(j), systMCclosureMax);
 548
 549       // FD (assume full correlations)
 550       errsyst=TMath::Abs(fhDplus->GetBinContent(j)*fSystDplus->GetHistoFDmin()->GetBinContent(j));
 551       systFDmin+=errsyst*weight;
 552       errsyst=fhDplus->GetBinContent(j)*fSystDplus->GetHistoFDmax()->GetBinContent(j);
 553       systFDmax+=errsyst*weight;
 554
 555       printf(" Dplus feeddown the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoFDmin()->GetBinContent(j), systMCclosureMin);
 556       printf(" Dplus feeddown the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoFDmax()->GetBinContent(j), systMCclosureMax);
 557
 558       printf("Dplus the value is: %f, weight: %f \n",value, weight);
 559
 560     }
 561
 562     value/=sumweights;
 563     errStatValue/=sumweights*sumweights;// was: ((Double_t)nmeson*sumweights);
 564     errSystYieldValue/=sumweights*sumweights;// was: ((Double_t)nmeson*sumweights);
 565     errSystBkgValue/=sumweights*sumweights; // was: ((Double_t)nmeson*sumweights);
 566     fhDaverage->SetBinContent(j,value);
 567     fhDaverage->SetBinError(j,TMath::Sqrt(errStatValue));
 568
 569     // Settting unc on yield and back sub
 570     fSystDaverage->GetHistoYieldUnc()->SetBinContent(j,TMath::Sqrt(errSystYieldValue)/TMath::Abs(value));
 571     fSystDaverage->GetHistoBackSubUncMin()->SetBinContent(j,-TMath::Sqrt(errSystBkgValue)/TMath::Abs(value));
 572     fSystDaverage->GetHistoBackSubUncMax()->SetBinContent(j,TMath::Sqrt(errSystBkgValue)/TMath::Abs(value));
 573
 574
 575
 576     // MCcorrections, associated tracks
 577     if(systMCcorrectionsMin>=0.){
 578       systMCcorrectionsMin=systMCcorrectionsMin/sumweights;
 579     }
 580     else {
 581       Printf("Negative sum in error calculation \n");
 582       return;
 583     }
 584
 585     if(systMCcorrectionsMax>=0.){
 586       systMCcorrectionsMax=systMCcorrectionsMax/sumweights;
 587     }
 588     else {
 589       Printf("Negative sum in error calculation \n");
 590       return;
 591     }
 592
 593
 594     // MCDefficiency
 595     if(systMCDefficiencyMin>=0.){
 596       systMCDefficiencyMin=TMath::Sqrt(systMCDefficiencyMin)/sumweights;
 597     }
 598     else {
 599       Printf("Negative sum in error calculation \n");
 600       return;
 601     }
 602
 603     if(systMCDefficiencyMax>=0.){
 604       systMCDefficiencyMax=TMath::Sqrt(systMCDefficiencyMax)/sumweights;
 605     }
 606     else {
 607       Printf("Negative sum in error calculation \n");
 608       return;
 609     }
 610
 611     // SecContamination
 612     if(systSecContaminationMin>=0.){
 613       systSecContaminationMin=systSecContaminationMin/sumweights;
 614     }
 615     else {
 616       Printf("Negative sum in error calculation \n");
 617       return;
 618     }
 619
 620     if(systSecContaminationMax>=0.){
 621       systSecContaminationMax=systSecContaminationMax/sumweights;
 622     }
 623     else {
 624       Printf("Negative sum in error calculation \n");
 625       return;
 626     }
 627
 628     // MCclosureTest
 629     if(systMCclosureMin>=0.){
 630       systMCclosureMin=systMCclosureMin/sumweights;
 631     }
 632     else {
 633       Printf("Negative sum in error calculation \n");
 634       return;
 635     }
 636
 637     if(systMCclosureMax>=0.){
 638       systMCclosureMax=systMCclosureMax/sumweights;
 639     }
 640     else {
 641       Printf("Negative sum in error calculation \n");
 642       return;
 643     }
 644
 645
 646
 647     // Feed down
 648     if(systFDmin>=0.){
 649       systFDmin=systFDmin/sumweights;
 650     }
 651     else {
 652       Printf("Negative sum in error calculation \n");
 653       return;
 654     }
 655
 656     if(systFDmax>=0.){
 657       systFDmax=systFDmax/sumweights;
 658     }
 659     else {
 660       Printf("Negative sum in error calculation \n");
 661       return;
 662     }
 663
 664     value=TMath::Abs(value);// protection to avoid flipping of min syst histo sign
 665
 666     fSystDaverage->GetHistoMCcorrectionsMin()->SetBinContent(j,-systMCcorrectionsMin/value);
 667     fSystDaverage->GetHistoMCcorrectionsMax()->SetBinContent(j,systMCcorrectionsMax/value);
 668
 669     fSystDaverage->GetHistoMCDefficiencyMin()->SetBinContent(j,-systMCDefficiencyMin/value);
 670     fSystDaverage->GetHistoMCDefficiencyMax()->SetBinContent(j,systMCDefficiencyMax/value);
 671
 672     fSystDaverage->GetHistoSecContaminationMin()->SetBinContent(j,-systSecContaminationMin/value);
 673     fSystDaverage->GetHistoSecContaminationMax()->SetBinContent(j,systSecContaminationMax/value);
 674
 675     fSystDaverage->GetHistoMCclosureTestMin()->SetBinContent(j,-systMCclosureMin/value);
 676     fSystDaverage->GetHistoMCclosureTestMax()->SetBinContent(j,systMCclosureMax/value);
 677
 678     fSystDaverage->GetHistoFDmin()->SetBinContent(j,-systFDmin/value);
 679     fSystDaverage->GetHistoFDmax()->SetBinContent(j,systFDmax/value);
 680
 681   }
 682   fSystDaverage->BuildTotalUncHisto();
 683   fSystDaverage->BuildTotalNonFDUncHisto();
 684   fSystDaverage->BuildGraphsUnc(fhDaverage);
 685   fSystDaverage->BuildGraphsRelUnc();
 686
 687   return ;
 688
 689 }
 690
 691
 692 TH1D* AliHFDmesonCorrAverage::ReflectHisto(TH1D *h){
 693
 694   TH1D *h2=new TH1D(Form("%sReflected",h->GetName()),Form("%sReflected",h->GetName()),h->GetNbinsX()/2.,0.,TMath::Pi());
 695   for(Int_t j=1;j<=h->GetNbinsX();j++){
 696     Double_t x=h->GetBinCenter(j);
 697     Double_t y0=h->GetBinContent(j);
 698     Double_t ey0=h->GetBinError(j);
 699     Int_t j2;
 700     if(x>0&&x<TMath::Pi())j2=h2->FindBin(x);
 701     else if(x<0)j2=h2->FindBin(-1.*x);
 702     else if(x>TMath::Pi())j2=h2->FindBin(2.*TMath::Pi()-x);
 703     else {
 704       printf("Point %d excluded \n",j);
 705       continue;
 706     }
 707     Double_t y=h2->GetBinContent(j2);
 708     Double_t ey=h2->GetBinError(j2);
 709     h2->SetBinContent(j2,y+y0);
 710     h2->SetBinError(j2,TMath::Sqrt(ey0*ey0+ey*ey));
 711   }
 712
 713   return h2;
 714 }
 715
 716 void AliHFDmesonCorrAverage::SetWeights(){
 717
 718
 719   if(fincludeDzero){
 720     TH1D *hYieldUnc=fSystDzero->GetHistoYieldUnc();
 721     TH1D *hBkgUncMin=fSystDzero->GetHistoBackSubUncMin();
 722     TH1D *hBkgUncMax=fSystDzero->GetHistoBackSubUncMax();
 723     if(fweightsDzeroStat)delete fweightsDzeroStat;
 724     if(fweightsDzeroSystYield)delete fweightsDzeroSystYield;
 725     if(fweightsDzeroSystBkg)delete fweightsDzeroSystBkg;
 726     fweightsDzeroStat=new Double_t[fnbinsphi];
 727     fweightsDzeroSystYield=new Double_t[fnbinsphi];
 728     fweightsDzeroSystBkg=new Double_t[fnbinsphi];
 729     //    fhGlobalWeightDzero=new TH1F("fhGlobalWeightDzero","fhGlobalWeightDzero",fnbinsphi
 730     for(Int_t j=0;j<fnbinsphi;j++){
 731       if(fmethod==10){
 732         fweightsDzeroStat[j]=1./(fhDzero->GetBinError(j+1)*fhDzero->GetBinError(j+1));
 733         fweightsDzeroSystYield[j]=1./(hYieldUnc->GetBinContent(j+1)*fhDzero->GetBinContent(j+1)*hYieldUnc->GetBinContent(j+1)*fhDzero->GetBinContent(j+1));
 734         //for asymmetric uncertainties...
 735         if(TMath::Abs(hBkgUncMin->GetBinContent(j+1)) > TMath::Abs(hBkgUncMax->GetBinContent(j+1))) fweightsDzeroSystBkg[j]=1./(hBkgUncMin->GetBinContent(j+1)*fhDzero->GetBinContent(j+1)*hBkgUncMin->GetBinContent(j+1)*fhDzero->GetBinContent(j+1));
 736         else fweightsDzeroSystBkg[j]=1./(hBkgUncMax->GetBinContent(j+1)*fhDzero->GetBinContent(j+1)*hBkgUncMax->GetBinContent(j+1)*fhDzero->GetBinContent(j+1));
 737       }
 738       else{
 739         Printf("This option for the calculation of the average has not been implemented yet");
 740         return;
 741       }
 742     }
 743   }
 744
 745   if(fincludeDstar){
 746     TH1D *hYieldUnc=fSystDstar->GetHistoYieldUnc();
 747     TH1D *hBkgUncMin=fSystDstar->GetHistoBackSubUncMin();
 748     TH1D *hBkgUncMax=fSystDstar->GetHistoBackSubUncMax();
 749     if(fweightsDstarStat)delete fweightsDstarStat;
 750     if(fweightsDstarSystYield)delete fweightsDstarSystYield;
 751     if(fweightsDstarSystBkg)delete fweightsDstarSystBkg;
 752     fweightsDstarStat=new Double_t[fnbinsphi];
 753     fweightsDstarSystYield=new Double_t[fnbinsphi];
 754     fweightsDstarSystBkg=new Double_t[fnbinsphi];
 755     for(Int_t j=0;j<fnbinsphi;j++){
 756       if(fmethod==10){
 757         fweightsDstarStat[j]=1./(fhDstar->GetBinError(j+1)*fhDstar->GetBinError(j+1));
 758         fweightsDstarSystYield[j]=1./(hYieldUnc->GetBinContent(j+1)*fhDstar->GetBinContent(j+1)*hYieldUnc->GetBinContent(j+1)*fhDstar->GetBinContent(j+1));
 759         //for asymmetric uncertainties...
 760         if(TMath::Abs(hBkgUncMin->GetBinContent(j+1)) > TMath::Abs(hBkgUncMax->GetBinContent(j+1))) fweightsDstarSystBkg[j]=1./(hBkgUncMin->GetBinContent(j+1)*fhDstar->GetBinContent(j+1)*hBkgUncMin->GetBinContent(j+1)*fhDstar->GetBinContent(j+1));
 761         else fweightsDstarSystBkg[j]=1./(hBkgUncMax->GetBinContent(j+1)*fhDstar->GetBinContent(j+1)*hBkgUncMax->GetBinContent(j+1)*fhDstar->GetBinContent(j+1));
 762       }
 763       else{
 764         Printf("This option for the calculation of the average has not been implemented yet");
 765         return;
 766       }
 767     }
 768   }
 769
 770   if(fincludeDplus){
 771     TH1D *hYieldUnc=fSystDplus->GetHistoYieldUnc();
 772     TH1D *hBkgUncMin=fSystDplus->GetHistoBackSubUncMin();
 773     TH1D *hBkgUncMax=fSystDplus->GetHistoBackSubUncMax();
 774     if(fweightsDplusStat)delete fweightsDplusStat;
 775     if(fweightsDplusSystYield)delete fweightsDplusSystYield;
 776     if(fweightsDplusSystBkg)delete fweightsDplusSystBkg;
 777     fweightsDplusStat=new Double_t[fnbinsphi];
 778     fweightsDplusSystYield=new Double_t[fnbinsphi];
 779     fweightsDplusSystBkg=new Double_t[fnbinsphi];
 780     for(Int_t j=0;j<fnbinsphi;j++){
 781       if(fmethod==10){
 782         fweightsDplusStat[j]=1./(fhDplus->GetBinError(j+1)*fhDplus->GetBinError(j+1));
 783         fweightsDplusSystYield[j]=1./(hYieldUnc->GetBinContent(j+1)*fhDplus->GetBinContent(j+1)*hYieldUnc->GetBinContent(j+1)*fhDplus->GetBinContent(j+1));
 784         //for asymmetric uncertainties...
 785         if(TMath::Abs(hBkgUncMin->GetBinContent(j+1)) > TMath::Abs(hBkgUncMax->GetBinContent(j+1))) fweightsDplusSystBkg[j]=1./(hBkgUncMin->GetBinContent(j+1)*fhDplus->GetBinContent(j+1)*hBkgUncMin->GetBinContent(j+1)*fhDplus->GetBinContent(j+1));
 786         else fweightsDplusSystBkg[j]=1./(hBkgUncMax->GetBinContent(j+1)*fhDplus->GetBinContent(j+1)*hBkgUncMax->GetBinContent(j+1)*fhDplus->GetBinContent(j+1));
 787       }
 788       else{
 789         Printf("This option for the calculation of the average has not been implemented yet");
 790         return;
 791       }
 792     }
 793   }
 794
 795
 796
 797
 798 }
 799
 800