PWG/muondep/AliAnalysisMuMuResult.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 // $Id$
  18
  19 #include "AliAnalysisMuMuResult.h"
  20
  21 ClassImp(AliAnalysisMuMuResult)
  22
  23 #include "TF1.h"
  24 #include "TFitResult.h"
  25 #include "TH1.h"
  26 #include "TH2.h"
  27 #include "THashList.h"
  28 #include "TLine.h"
  29 #include "TList.h"
  30 #include "TMap.h"
  31 #include "TMath.h"
  32 #include "TObjArray.h"
  33 #include "TParameter.h"
  34 #include "AliAnalysisMuMuBinning.h"
  35 #include "AliHistogramCollection.h"
  36 #include "AliLog.h"
  37 #include <map>
  38
  39 const std::map<std::string,Double_t>& MassMap()
  40 {
  41   static std::map<std::string,Double_t> massMap;
  42   // not super elegant, but that way we do not depend on TDatabasePDG and thus
  43   // can decide on our particle naming
  44   if (massMap.empty())
  45   {
  46     massMap["Jpsi"] = 3.096916e+00;
  47     massMap["PsiPrime"] = 3.68609e+00;
  48     massMap["Upsilon"] = 9.46030e+00;
  49     massMap["UpsilonPrime"] = 1.00233e+01;
  50   }
  51   return massMap;
  52 }
  53
  54
  55 Double_t ErrorAB(Double_t a, Double_t aerr, Double_t b, Double_t berr)
  56 {
  57   Double_t e(0.0);
  58
  59   if ( TMath::Abs(a) > 1E-12 )
  60   {
  61     e += (aerr*aerr)/(a*a);
  62   }
  63
  64   if ( TMath::Abs(b) > 1E-12 )
  65   {
  66     e += (berr*berr)/(b*b);
  67   }
  68
  69   return TMath::Sqrt(e);
  70 }
  71
  72 Double_t funcCB(Double_t* xx, Double_t* par)
  73 {
  74   Double_t N = par[0];
  75   Double_t alpha = par[1];
  76   Double_t n = par[2];
  77   Double_t mean = par[3];
  78   Double_t sigma = par[4];
  79
  80   Double_t x = xx[0];
  81
  82   Double_t A = TMath::Power(n/TMath::Abs(alpha),n)*TMath::Exp(-0.5*alpha*alpha);
  83   Double_t B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
  84
  85   Double_t y = ( TMath::Abs(sigma) > 1E-12 ? (x-mean)/sigma : 0 );
  86
  87   if ( y > alpha*-1.0 )
  88   {
  89     return N*TMath::Exp(-0.5*y*y);
  90   }
  91   else
  92   {
  93     return N*A*TMath::Power(B-y,-n);
  94   }
  95 }
  96
  97 Double_t funcJpsiGCBE(Double_t* xx, Double_t* par)
  98 {
  99   Double_t x = xx[0];
 100
 101   Double_t g = par[0]*TMath::Gaus(x,par[1],par[2]);
 102
 103   Double_t jpsi = funcCB(xx,par+3);
 104
 105   Double_t expo = par[8]*TMath::Exp(par[9]*x);
 106
 107   return g+expo+jpsi;
 108 }
 109
 110 Double_t funcJpsiPsiPrimeCustom(Double_t* xx, Double_t* par)
 111 {
 112   Double_t N = par[0];
 113   Double_t alpha = par[1];
 114   Double_t n = par[2];
 115   Double_t mean = par[3];
 116   Double_t sigma = par[4];
 117   Double_t alphaprime = par[5];
 118   Double_t nprime = par[6];
 119
 120   Double_t x = xx[0];
 121
 122   Double_t A = TMath::Power(n/TMath::Abs(alpha),n)*TMath::Exp(-0.5*alpha*alpha);
 123   Double_t B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
 124   Double_t C = TMath::Power(nprime/TMath::Abs(alphaprime),nprime)*TMath::Exp(-0.5*alphaprime*alphaprime);
 125   Double_t D = nprime/TMath::Abs(alphaprime) - TMath::Abs(alphaprime);
 126
 127   Double_t y = ( TMath::Abs(sigma) > 1E-12 ? (x-mean)/sigma : 0 );
 128
 129   Double_t cb(0);
 130
 131   if ( y > alphaprime )
 132   {
 133     cb = N*C*TMath::Power(D+y,-nprime);
 134   }
 135   else if ( y > alpha*-1.0 )
 136   {
 137     cb = N*TMath::Exp(-0.5*y*y);
 138   }
 139   else
 140   {
 141     cb = N*A*TMath::Power(B-y,-n);
 142   }
 143
 144   if ( x < mean )
 145   {
 146     return cb + par[7] + par[8]*x; // gaus + pol1
 147   }
 148   else
 149   {
 150     Double_t yprime = (x-par[10])/par[11];
 151     return cb + par[9]*TMath::Exp(-0.5*yprime*yprime)+par[12]*TMath::Exp(-par[13]*x);
 152     // gaus (j/psi) + gaus (psi') + expo
 153   }
 154 }
 155
 156
 157 Double_t funcCB2(Double_t* xx, Double_t* par)
 158 {
 159   Double_t N = par[0];
 160   Double_t alpha = par[1];
 161   Double_t n = par[2];
 162   Double_t mean = par[3];
 163   Double_t sigma = par[4];
 164   Double_t alphaprime = par[5];
 165   Double_t nprime = par[6];
 166
 167   Double_t x = xx[0];
 168
 169   Double_t A = TMath::Power(n/TMath::Abs(alpha),n)*TMath::Exp(-0.5*alpha*alpha);
 170   Double_t B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
 171   Double_t C = TMath::Power(nprime/TMath::Abs(alphaprime),nprime)*TMath::Exp(-0.5*alphaprime*alphaprime);
 172   Double_t D = nprime/TMath::Abs(alphaprime) - TMath::Abs(alphaprime);
 173
 174   Double_t y = ( TMath::Abs(sigma) > 1E-12 ? (x-mean)/sigma : 0 );
 175
 176   if ( y > alphaprime )
 177   {
 178     return N*C*TMath::Power(D+y,-nprime);
 179   }
 180   else if ( y > alpha*-1.0 )
 181   {
 182     return N*TMath::Exp(-0.5*y*y);
 183   }
 184   else
 185   {
 186     return N*A*TMath::Power(B-y,-n);
 187   }
 188 }
 189
 190 Double_t funcJpsiPsiPrime(Double_t* xx, Double_t* par)
 191 {
 192   Double_t jpsi = funcCB(xx,par);
 193   Double_t psiprime = funcCB2(xx,par+5);
 194
 195   int n = 10;
 196   Double_t x = xx[0];
 197
 198   Double_t e1 = par[n]*TMath::Exp(par[n+1]*x);
 199   Double_t e2 = par[n+2]*TMath::Exp(par[n+3]*x);
 200
 201   Double_t e = e1;
 202
 203   if ( x > par[3] ) e=e2;
 204
 205   return jpsi+psiprime+e;
 206 }
 207
 208 Double_t funcJpsiCBE(Double_t* xx, Double_t* par)
 209 {
 210   // CB + expo
 211
 212   Double_t jpsi = funcCB(xx,par);
 213
 214   Double_t x = xx[0];
 215
 216   Double_t e1 = par[5]*TMath::Exp(par[6]*x);
 217
 218   return jpsi+e1;
 219 }
 220
 221
 222 Double_t funcJpsiPCBE(Double_t* xx, Double_t* par)
 223 {
 224   Double_t x = xx[0];
 225
 226   Double_t pol2 = par[0] + par[1]*x + par[2]*x*x;
 227
 228   Double_t jpsi = funcCB(xx,par+3);
 229
 230   Double_t expo = par[8]*TMath::Exp(par[9]*x);
 231
 232   return pol2+jpsi+expo;
 233 }
 234
 235 Double_t funcJpsiECBE(Double_t* xx, Double_t* par)
 236 {
 237   // CB + expo
 238
 239   Double_t jpsi = funcCB(xx,par+2);
 240
 241   Double_t x = xx[0];
 242
 243   Double_t e1 = par[0]*TMath::Exp(par[1]*x);
 244
 245   Double_t e2 = par[7]*TMath::Exp(par[8]*x);
 246
 247   return e1+e2+jpsi;
 248 }
 249
 250 const char* NormalizeName(const char* name, const char* suffix)
 251 {
 252   TString str(Form("%s_%s",name,suffix));
 253
 254   str.ReplaceAll("-","_");
 255   str.ReplaceAll("/","%");
 256
 257   return str.Data();
 258 }
 259
 260 //------------------------------------------------------------------------------
 261 Double_t BackgroundVWG(Double_t *x, Double_t *par)
 262 {
 263   // gaussian variable width
 264   Double_t sigma = par[2]+par[3]*((x[0]-par[1])/par[1]);
 265   return par[0]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/(2.*sigma*sigma));
 266
 267 }
 268
 269 //------------------------------------------------------------------------------
 270 Double_t CrystalBallExtended(Double_t *x,Double_t *par)
 271 {
 272   //par[0] = Normalization
 273   //par[1] = mean
 274   //par[2] = sigma
 275   //par[3] = alpha
 276   //par[4] = n
 277   //par[5] = alpha'
 278   //par[6] = n'
 279
 280   Double_t t = (x[0]-par[1])/par[2];
 281   if (par[3] < 0) t = -t;
 282
 283   Double_t absAlpha = fabs((Double_t)par[3]);
 284   Double_t absAlpha2 = fabs((Double_t)par[5]);
 285
 286   if (t >= -absAlpha && t < absAlpha2) // gaussian core
 287   {
 288     return par[0]*(exp(-0.5*t*t));
 289   }
 290
 291   if (t < -absAlpha) //left tail
 292   {
 293     Double_t a =  TMath::Power(par[4]/absAlpha,par[4])*exp(-0.5*absAlpha*absAlpha);
 294     Double_t b = par[4]/absAlpha - absAlpha;
 295     return par[0]*(a/TMath::Power(b - t, par[4]));
 296   }
 297
 298   if (t >= absAlpha2) //right tail
 299   {
 300
 301     Double_t c =  TMath::Power(par[6]/absAlpha2,par[6])*exp(-0.5*absAlpha2*absAlpha2);
 302     Double_t d = par[6]/absAlpha2 - absAlpha2;
 303     return par[0]*(c/TMath::Power(d + t, par[6]));
 304   }
 305
 306   return 0. ;
 307 }
 308
 309 //------------------------------------------------------------------------------
 310 Double_t Gaus(Double_t *x, Double_t *par)
 311 {
 312   // gaussian
 313   return par[0]/TMath::Sqrt(2.*TMath::Pi())/par[2]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/(2.*par[2]*par[2]));
 314
 315 }
 316
 317 //------------------------------------------------------------------------------
 318 Double_t Exp(Double_t *x, Double_t *par)
 319 {
 320   // exponential
 321   return par[0]*TMath::Exp(par[1]*x[0]);
 322
 323 }
 324
 325 //------------------------------------------------------------------------------
 326 Double_t Pow(Double_t *x, Double_t *par)
 327 {
 328   // power law
 329   return par[0]*TMath::Power(x[0],par[1]);
 330
 331 }
 332
 333 //------------------------------------------------------------------------------
 334 Double_t fitFunctionVWG(Double_t *x, Double_t *par)
 335 {
 336   if (x[0] > 2.9 && x[0] < 3.3) TF1::RejectPoint();
 337   return BackgroundVWG(x, par);
 338 }
 339
 340 //------------------------------------------------------------------------------
 341 Double_t fitFunctionCB2VWG(Double_t *x, Double_t *par)
 342 {
 343   return BackgroundVWG(x, par) + CrystalBallExtended(x, &par[4]);
 344 }
 345
 346 //------------------------------------------------------------------------------
 347 Double_t func2CB2VWG(Double_t *x, Double_t *par)
 348 {
 349   Double_t par2[7] = {
 350     par[11],
 351     3.68609+(par[5]-3.096916)/3.096916*3.68609,
 352     par[6]/3.096916*3.68609,
 353     par[7],
 354     par[8],
 355     par[9],
 356     par[10]
 357   };
 358   return BackgroundVWG(x, par) + CrystalBallExtended(x, &par[4]) + CrystalBallExtended(x, par2);
 359 }
 360
 361 //_____________________________________________________________________________
 362 //_____________________________________________________________________________
 363 //_____________________________________________________________________________
 364 //_____________________________________________________________________________
 365 //_____________________________________________________________________________
 366
 367 //_____________________________________________________________________________
 368 AliAnalysisMuMuResult::AliAnalysisMuMuResult(TRootIOCtor* /*io*/) :
 369 TNamed("",""),
 370 fNofRuns(),
 371 fNofTriggers(-1),
 372 fMinv(0x0),
 373 fBin(),
 374 fSubResults(0x0),
 375 fMap(0x0),
 376 fMother(0x0),
 377 fKeys(0x0),
 378 fWeight(0.0)
 379 {
 380 }
 381
 382 //_____________________________________________________________________________
 383 AliAnalysisMuMuResult::AliAnalysisMuMuResult(const TH1& hminv)
 384 :
 385   TNamed("",""),
 386   fNofRuns(1),
 387   fNofTriggers(-1),
 388   fMinv(0x0),
 389   fBin(),
 390   fSubResults(0x0),
 391   fMap(0x0),
 392   fMother(0x0),
 393 fKeys(0x0),
 394 fWeight(0.0)
 395
 396 {
 397   SetMinv(hminv);
 398 }
 399
 400 //_____________________________________________________________________________
 401 AliAnalysisMuMuResult::AliAnalysisMuMuResult(const TH1& hminv,
 402                                              const char* fitType,
 403                                              Int_t nrebin)
 404 :
 405 TNamed(Form("%s-%d",fitType,nrebin),""),
 406 fNofRuns(1),
 407 fNofTriggers(-1),
 408 fMinv(0x0),
 409 fBin(),
 410 fSubResults(0x0),
 411 fMap(0x0),
 412 fMother(0x0),
 413 fKeys(0x0),
 414 fWeight(0.0)
 415
 416 {
 417   SetMinv(hminv);
 418 }
 419
 420 //_____________________________________________________________________________
 421 AliAnalysisMuMuResult::AliAnalysisMuMuResult(const TH1& hminv,
 422                                              const char* triggerName,
 423                                              const char* eventSelection,
 424                                              const char* pairSelection,
 425                                              const char* centSelection,
 426                                              const AliAnalysisMuMuBinning::Range& bin)
 427 :
 428 TNamed(Form("%s-%s-%s-%s",triggerName,eventSelection,pairSelection,centSelection),""),
 429 fNofRuns(1),
 430 fNofTriggers(-1),
 431 fMinv(0x0),
 432 fBin(bin),
 433 fSubResults(0x0),
 434 fMap(0x0),
 435 fMother(0x0),
 436 fKeys(0x0),
 437 fWeight(0.0)
 438
 439 {
 440   SetMinv(hminv);
 441 }
 442
 443 //_____________________________________________________________________________
 444 AliAnalysisMuMuResult::AliAnalysisMuMuResult(const AliAnalysisMuMuResult& rhs)
 445 :
 446 TNamed(rhs),
 447 fNofRuns(rhs.NofRuns()),
 448 fNofTriggers(rhs.NofTriggers()),
 449 fMinv(0x0),
 450 fBin(rhs.Bin()),
 451 fSubResults(0x0),
 452 fMap(0x0),
 453 fMother(0x0),
 454 fKeys(0x0),
 455 fWeight(rhs.fWeight)
 456 {
 457   /// copy ctor
 458   /// Note that the mother is lost
 459   /// fKeys remains 0x0 so it will be recomputed if need be
 460
 461   if ( rhs.fMinv )
 462   {
 463     fMinv = static_cast<TH1*>(rhs.fMinv->Clone());
 464   }
 465
 466   if (rhs.fSubResults)
 467   {
 468     fSubResults = static_cast<TObjArray*>(rhs.fSubResults->Clone());
 469   }
 470
 471   if ( rhs.fMap )
 472   {
 473     fMap = static_cast<TMap*>(rhs.fMap->Clone());
 474   }
 475
 476 }
 477
 478 //_____________________________________________________________________________
 479 AliAnalysisMuMuResult& AliAnalysisMuMuResult::operator=(const AliAnalysisMuMuResult& rhs)
 480 {
 481   if (this!=&rhs)
 482   {
 483     delete fMinv;
 484     delete fMap;
 485     delete fSubResults;
 486
 487     fMinv = 0x0;
 488     fMap = 0x0;
 489     fSubResults = 0x0;
 490     fKeys = 0x0;
 491
 492     if ( rhs.fMinv )
 493     {
 494       fMinv = static_cast<TH1*>(rhs.fMinv->Clone());
 495     }
 496
 497     if (rhs.fSubResults)
 498     {
 499       fSubResults = static_cast<TObjArray*>(rhs.fSubResults->Clone());
 500     }
 501
 502     if ( rhs.fMap )
 503     {
 504       fMap = static_cast<TMap*>(rhs.fMap->Clone());
 505     }
 506
 507     static_cast<TNamed&>(*this)=rhs;
 508
 509     fNofRuns = rhs.NofRuns();
 510     fNofTriggers = rhs.NofTriggers();
 511     fBin = rhs.Bin();
 512     fWeight = rhs.fWeight;
 513   }
 514
 515   return *this;
 516 }
 517
 518 //_____________________________________________________________________________
 519 AliAnalysisMuMuResult::~AliAnalysisMuMuResult()
 520 {
 521   delete fMap;
 522   delete fMinv;
 523   delete fSubResults;
 524   delete fKeys;
 525 }
 526
 527 //_____________________________________________________________________________
 528 const AliAnalysisMuMuBinning::Range& AliAnalysisMuMuResult::Bin() const
 529 {
 530   if ( !Mother() ) return fBin;
 531   else return Mother()->Bin();
 532 }
 533
 534 //_____________________________________________________________________________
 535 TObject* AliAnalysisMuMuResult::Clone(const char* /*newname*/) const
 536 {
 537   return new AliAnalysisMuMuResult(*this);
 538 }
 539
 540 //_____________________________________________________________________________
 541 Bool_t AliAnalysisMuMuResult::Correct(const AliAnalysisMuMuResult& other, const char* particle, const char* subResultName)
 542 {
 543   /// Assuming other has an AccxEff entry, correct this value by AccxEff of other
 544
 545   if ( HasValue(Form("Nof%s",particle)) )
 546   {
 547     if (!other.HasValue(Form("AccEff%s",particle),subResultName))
 548     {
 549       AliError(Form("Cannot correct as I do not find the AccEff%s value (subResultName=%s)!",particle,subResultName));
 550       return kFALSE;
 551     }
 552
 553     Double_t acc = other.GetValue(Form("AccEff%s",particle),subResultName);
 554     Double_t value = 0.0;
 555
 556     if ( acc > 0 ) value = GetValue(Form("Nof%s",particle)) / acc;
 557
 558     Double_t error = ErrorAB( GetValue(Form("Nof%s",particle)),
 559                               GetErrorStat(Form("Nof%s",particle)),
 560                               other.GetValue(Form("AccEff%s",particle),subResultName),
 561                               other.GetErrorStat(Form("AccEff%s",particle),subResultName) );
 562
 563     Set(Form("CorrNof%s",particle),value,error*value);
 564
 565     return kTRUE;
 566   }
 567
 568   return kFALSE;
 569 }
 570
 571 //_____________________________________________________________________________
 572 Double_t AliAnalysisMuMuResult::CountParticle(const TH1& hminv, const char* particle, Double_t sigma)
 573 {
 574   /// Count the number of entries in an invariant mass range
 575
 576   const std::map<std::string,Double_t>& m = MassMap();
 577
 578   std::map<std::string,Double_t>::const_iterator it = m.find(particle);
 579
 580   if ( it == m.end() )
 581   {
 582     AliErrorClass(Form("Don't know the mass of particle %s",particle));
 583     return 0.0;
 584   }
 585
 586   Double_t mass = it->second;
 587
 588   if ( sigma < 0 )
 589   {
 590     AliDebugClass(1,Form("Oups. Got a sigma of %e for particle %s !",sigma,particle));
 591     return hminv.Integral();
 592   }
 593
 594   TAxis* x = hminv.GetXaxis();
 595
 596   Int_t b1 = x->FindBin(mass-sigma);
 597   Int_t b2 = x->FindBin(mass+sigma);
 598
 599   AliDebugClass(1,Form("hminv getentries %e integral %e",hminv.GetEntries(),hminv.Integral(b1,b2)));
 600
 601   return hminv.Integral(b1,b2);
 602 }
 603
 604 /*
 605 //_____________________________________________________________________________
 606 void AliAnalysisMuMuResult::FitJpsiPsiPrimeCustom(TH1& h)
 607 {
 608   std::cout << "Fit with jpsi + psiprime (custom)" << std::endl;
 609
 610   const Double_t xmin(1.5);
 611   const Double_t xmax(8.0);
 612
 613   fitTotal = new TF1("fitTotal",funcJpsiPsiPrimeCustom,xmin,xmax,14);
 614   fitTotal->SetLineColor(4);
 615
 616   fitTotal->SetParName(0,"cstecb");
 617   fitTotal->SetParName(1,"alpha");
 618   fitTotal->SetParName(2,"n");
 619   fitTotal->SetParName(3,"meanjpsi");
 620   fitTotal->SetParName(4,"sigmajpsi");
 621   fitTotal->SetParName(5,"alphaprime");
 622   fitTotal->SetParName(6,"nprime");
 623   fitTotal->SetParName(7,"cstepol1");
 624   fitTotal->SetParName(8,"slopepol1");
 625   fitTotal->SetParName(9,"cstegaus");
 626   fitTotal->SetParName(10,"meanpsiprime");
 627   fitTotal->SetParName(11,"sigmapsiprime");
 628   fitTotal->SetParName(12,"csteexpo");
 629   fitTotal->SetParName(13,"slopeexpo");
 630
 631   fitTotal->SetParameter( 0,1);
 632
 633   const char* fitOption = "SQBR+";
 634   const Double_t alphaMC = 0.936;
 635   const Double_t nMC = 4.44;
 636   const Double_t alphaprimeMC = 1.60;
 637   const Double_t nprimeMC = 3.23;
 638
 639   TF1* fcb = new TF1("cb",funcCB2,2.9,3.3,7);
 640   fcb->SetParameters(1,1.0,4.0,3.1,0.1,1.5,3);
 641
 642   fcb->SetParLimits(3,3,4);
 643   fcb->SetParLimits(4,0,1);
 644
 645   fcb->FixParameter(1,alphaMC);
 646   fcb->FixParameter(2,nMC);
 647   fcb->FixParameter(5,alphaprimeMC);
 648   fcb->FixParameter(6,nprimeMC);
 649
 650   TFitResultPtr rcb = h.Fit(fcb,fitOption,"",2.9,3.3);
 651
 652   if (!rcb.Get())
 653   {
 654     return;
 655   }
 656
 657   fitTotal->SetParameter(0,rcb->Parameter(0));
 658   fitTotal->SetParameter(1,rcb->Parameter(1)); fitTotal->SetParLimits(1,0,10); // alpha
 659   fitTotal->SetParameter(2,rcb->Parameter(2)); fitTotal->SetParLimits(2,1,10); // n
 660   fitTotal->SetParameter(3,rcb->Parameter(3)); fitTotal->SetParLimits(3,3.0,3.5); // mean
 661   fitTotal->SetParameter(4,rcb->Parameter(4)); fitTotal->SetParLimits(4,0,1); // sigma
 662   fitTotal->SetParameter(5,rcb->Parameter(5)); fitTotal->SetParLimits(1,0,10); // alphaprime
 663   fitTotal->SetParameter(6,rcb->Parameter(6)); fitTotal->SetParLimits(2,1,10); // nprime
 664
 665   fitTotal->FixParameter(1,alphaMC);
 666   fitTotal->FixParameter(2,nMC);
 667   fitTotal->FixParameter(5,alphaprimeMC);
 668   fitTotal->FixParameter(6,nprimeMC);
 669
 670   TF1* fge = new TF1("fge","gaus(0)+expo(3)",3.5,4.4);
 671   fge->SetParameters(1,3.6,0.25,1,1);
 672   TFitResultPtr rpsiprime = h.Fit(fge,fitOption,"",3.5,4.4);
 673
 674   if (static_cast<int>(rpsiprime))
 675   {
 676     AliInfo("Will fix psiprime parameters");
 677     fitTotal->FixParameter(9,0);
 678     fitTotal->FixParameter(10,3.7);
 679     fitTotal->FixParameter(11,0.1);
 680   }
 681   else
 682   {
 683     fitTotal->SetParameter(10,rpsiprime->Parameter(1)); fitTotal->SetParLimits(10,3.5,3.8); // mean'
 684     fitTotal->SetParameter(11,rpsiprime->Parameter(2)); fitTotal->SetParLimits(11,0.05,0.7); // sigma'
 685   }
 686
 687   TFitResultPtr rpol1 = h.Fit("pol1",fitOption,"",1.5,2.5);
 688   fitTotal->SetParameter( 7,rpol1->Parameter(0));
 689   fitTotal->SetParameter( 8,rpol1->Parameter(1));
 690
 691   TFitResultPtr rexpo = h.Fit("expo",fitOption,"",4.5,7.0);
 692   fitTotal->SetParameter(12,rexpo->Parameter(0));
 693   fitTotal->SetParameter(13,rexpo->Parameter(1));
 694
 695
 696   TFitResultPtr r = h.Fit(fitTotal,fitOption,"",1.5,7);
 697
 698   TF1* signal = new TF1("signal","gaus",2,6);
 699   signal->SetParameters(fitTotal->GetParameter(0),
 700                         fitTotal->GetParameter(3),
 701                         fitTotal->GetParameter(4));
 702
 703   TF1* signalPrime = new TF1("signalPrime","gaus",2,6);
 704   signalPrime->SetParameters(fitTotal->GetParameter(9),
 705                              fitTotal->GetParameter(10),
 706                              fitTotal->GetParameter(11));
 707
 708   Double_t gausParameters[3];
 709   Double_t covarianceMatrix[3][3];
 710   Double_t gausParametersPrime[3];
 711   Double_t covarianceMatrixPrime[3][3];
 712
 713   covarianceMatrix[0][0] = (r->GetCovarianceMatrix())(0,0);
 714   covarianceMatrix[1][0] = (r->GetCovarianceMatrix())(3,0);
 715   covarianceMatrix[2][0] = (r->GetCovarianceMatrix())(4,0);
 716   covarianceMatrix[0][1] = (r->GetCovarianceMatrix())(0,3);
 717   covarianceMatrix[0][2] = (r->GetCovarianceMatrix())(0,4);
 718
 719   for ( int iy = 1; iy < 3; ++iy )
 720   {
 721     for ( int ix = 1; ix < 3; ++ix )
 722     {
 723       covarianceMatrix[ix][iy] = (r->GetCovarianceMatrix())(ix+2,iy+2);
 724     }
 725   }
 726
 727   gausParameters[0] = fitTotal->GetParameter(0);
 728   gausParameters[1] = fitTotal->GetParameter(3);
 729   gausParameters[2] = fitTotal->GetParameter(4);
 730
 731   gausParametersPrime[0] = fitTotal->GetParameter(9);
 732   gausParametersPrime[1] = fitTotal->GetParameter(10);
 733   gausParametersPrime[2] = fitTotal->GetParameter(11);
 734
 735   covarianceMatrixPrime[0][0] = (r->GetCovarianceMatrix())(9,9);
 736   covarianceMatrixPrime[1][0] = (r->GetCovarianceMatrix())(10,9);
 737   covarianceMatrixPrime[2][0] = (r->GetCovarianceMatrix())(11,9);
 738   covarianceMatrixPrime[0][1] = (r->GetCovarianceMatrix())(9,10);
 739   covarianceMatrixPrime[0][2] = (r->GetCovarianceMatrix())(9,11);
 740
 741   for ( int iy = 1; iy < 3; ++iy )
 742   {
 743     for ( int ix = 1; ix < 3; ++ix )
 744     {
 745       covarianceMatrixPrime[ix][iy] = (r->GetCovarianceMatrix())(ix+2,iy+2);
 746     }
 747   }
 748
 749   double n = signal->Integral(2,6)/h.GetBinWidth(10);
 750   double nerr = signal->IntegralError(2,6,&gausParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(10);
 751   Set("NofJpsi",n,nerr);
 752   Set("MeanJpsi",fitTotal->GetParameter(3),fitTotal->GetParError(3));
 753   Set("SigmaJpsi",fitTotal->GetParameter(4),fitTotal->GetParError(4));
 754
 755   double nprime = signalPrime->Integral(2,6)/h.GetBinWidth(10);
 756   double nerrprime = signalPrime->IntegralError(2,6,&gausParametersPrime[0],&covarianceMatrixPrime[0][0])/h.GetBinWidth(10);
 757   Set("NofPsiPrime",nprime,nerrprime);
 758   Set("MeanPsiPrime",fitTotal->GetParameter(10),fitTotal->GetParError(10));
 759   Set("SigmaPsiPrime",fitTotal->GetParameter(11),fitTotal->GetParError(11));
 760 }
 761 */
 762
 763 /*
 764 //_____________________________________________________________________________
 765 AliAnalysisMuMuResult::SubResult AliAnalysisMuMuResult::FitJpsiPsiPrimeCB(TH1& h)
 766 {
 767   std::cout << "Fit with jpsi + psiprime (CB) " << std::endl;
 768
 769   const Double_t xmin(1.5);
 770   const Double_t xmax(8.0);
 771
 772   fitTotal = new TF1("fitTotal",funcJpsiPsiPrime,xmin,xmax,14);
 773
 774 //  Double_t N = par[0];
 775 //  Double_t alpha = par[1];
 776 //  Double_t n = par[2];
 777 //  Double_t mean = par[3];
 778 //  Double_t sigma = par[4];
 779
 780   fitTotal->SetParameter( 0,1); // N
 781   fitTotal->FixParameter( 1,0.936); // alpha
 782   fitTotal->FixParameter( 2,4.44); // n
 783   fitTotal->SetParameter( 3,3.1); fitTotal->SetParLimits(3,3.0,3.2); // mean
 784   fitTotal->SetParameter( 4,0.07); fitTotal->SetParLimits(4,0.02,1); // sigma
 785
 786   fitTotal->SetParameter( 5,0.01); // N'
 787   fitTotal->FixParameter( 6,0.936); // alpha'
 788   fitTotal->FixParameter( 7,4.44); // n'
 789   fitTotal->SetParameter( 8,3.7); fitTotal->SetParLimits(8,3.5,3.8); // mean'
 790   fitTotal->SetParameter( 9,0.1); fitTotal->SetParLimits(9,0.02,1.0); // sigma'
 791
 792   fitTotal->SetParameter(10,h.GetMaximum());
 793   fitTotal->SetParameter(11,-1);
 794
 795   fitTotal->SetParameter(12,h.GetMaximum()/100);
 796   fitTotal->SetParameter(13,-1);
 797
 798   TFitResultPtr r = h.Fit(fitTotal,"SQBI","",1.5,6);
 799
 800 //  for ( int ix = 0; ix < fitTotal->GetNpar(); ++ix )
 801 //  {
 802 //    for ( int iy = 0; iy < fitTotal->GetNpar(); ++iy )
 803 //    {
 804 //      std::cout << Form("COV(%d,%d)=%e ",ix,iy,r->GetCovarianceMatrix()(ix,iy));
 805 //    }
 806 //    std::cout << std::endl;
 807 //  }
 808
 809
 810   TF1* signal = new TF1("signal","gaus",2,8);
 811
 812   signal->SetParameters(fitTotal->GetParameter(0),
 813                         fitTotal->GetParameter(3),
 814                         fitTotal->GetParameter(4));
 815
 816   TF1* signalPrime = new TF1("signalPrime","gaus",2,8);
 817
 818   signalPrime->SetParameters(fitTotal->GetParameter(0),
 819                              fitTotal->GetParameter(8),
 820                              fitTotal->GetParameter(9));
 821
 822   Double_t gausParameters[3];
 823   Double_t gausParametersPrime[3];
 824   Double_t covarianceMatrix[3][3];
 825   Double_t covarianceMatrixPrime[3][3];
 826
 827   gausParameters[0] = fitTotal->GetParameter(0);
 828   gausParameters[1] = fitTotal->GetParameter(3);
 829   gausParameters[2] = fitTotal->GetParameter(4);
 830
 831   covarianceMatrix[0][0] = (r->GetCovarianceMatrix())(0,0);
 832   covarianceMatrix[1][0] = (r->GetCovarianceMatrix())(3,0);
 833   covarianceMatrix[2][0] = (r->GetCovarianceMatrix())(4,0);
 834   covarianceMatrix[0][1] = (r->GetCovarianceMatrix())(0,3);
 835   covarianceMatrix[0][2] = (r->GetCovarianceMatrix())(0,4);
 836
 837   for ( int iy = 1; iy < 3; ++iy )
 838   {
 839     for ( int ix = 1; ix < 3; ++ix )
 840     {
 841       covarianceMatrix[ix][iy] = (r->GetCovarianceMatrix())(ix+2,iy+2);
 842     }
 843   }
 844
 845   gausParametersPrime[0] = fitTotal->GetParameter(0);
 846   gausParametersPrime[1] = fitTotal->GetParameter(8);
 847   gausParametersPrime[2] = fitTotal->GetParameter(9);
 848
 849   covarianceMatrixPrime[0][0] = (r->GetCovarianceMatrix())(0,0);
 850   covarianceMatrixPrime[1][0] = (r->GetCovarianceMatrix())(8,0);
 851   covarianceMatrixPrime[2][0] = (r->GetCovarianceMatrix())(9,0);
 852   covarianceMatrixPrime[0][1] = (r->GetCovarianceMatrix())(0,8);
 853   covarianceMatrixPrime[0][2] = (r->GetCovarianceMatrix())(0,9);
 854
 855   for ( int iy = 1; iy < 3; ++iy )
 856   {
 857     for ( int ix = 1; ix < 3; ++ix )
 858     {
 859       covarianceMatrixPrime[ix][iy] = (r->GetCovarianceMatrix())(ix+2,iy+2);
 860     }
 861   }
 862
 863   double n = signal->Integral(2,6)/h.GetBinWidth(10);
 864   double nerr = signal->IntegralError(2,6,&gausParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(10);
 865
 866   Set("NofJpsi",n,nerr);
 867   Set("MeanJpsi",fitTotal->GetParameter(3),fitTotal->GetParError(3));
 868   Set("SigmaJpsi",fitTotal->GetParameter(4),fitTotal->GetParError(4));
 869
 870   double nprime = signalPrime->Integral(2,6)/h.GetBinWidth(10);
 871   double nerrprime = signalPrime->IntegralError(2,6,&gausParametersPrime[0],&covarianceMatrixPrime[0][0])/h.GetBinWidth(10);
 872
 873   Set("NofPsiPrime",nprime,nerrprime);
 874   Set("MeanPsiPrime",fitTotal->GetParameter(8),fitTotal->GetParError(8));
 875   Set("SigmaPsiPrime",fitTotal->GetParameter(9),fitTotal->GetParError(9));
 876
 877 }
 878   */
 879
 880 //_____________________________________________________________________________
 881 AliAnalysisMuMuResult* AliAnalysisMuMuResult::FitJpsiGCBE(TH1& h)
 882 {
 883   std::cout << "Fit with jpsi alone (gaus + CB + expo)" << std::endl;
 884
 885   Int_t nrebin = fMinv->GetXaxis()->GetNbins() / h.GetXaxis()->GetNbins();
 886
 887   AliAnalysisMuMuResult* r = new AliAnalysisMuMuResult(h,"JPSIGCBE",nrebin);
 888
 889   TH1* hfit = r->Minv();
 890
 891   const Double_t xmin(1.0);
 892   const Double_t xmax(8.0);
 893
 894   TF1* fitTotal = new TF1("fitTotal",funcJpsiGCBE,xmin,xmax,10);
 895   fitTotal->SetParNames("cste","x0","sigma0","N","alpha","n","mean","sigma","expocste","exposlope");
 896
 897   fitTotal->SetParLimits(3,0,h.GetMaximum()*2); // N
 898
 899   const Double_t cbalpha(0.98);
 900   const Double_t cbn(5.2);
 901
 902   fitTotal->FixParameter(4,cbalpha);
 903   fitTotal->FixParameter(5,cbn);
 904
 905   fitTotal->SetParLimits(6,2.8,3.2); // mean
 906   fitTotal->SetParLimits(7,0.02,0.3); // sigma
 907
 908   TF1* fg = new TF1("fg","gaus",xmin,xmax);
 909
 910   hfit->Fit(fg,"","",0.75,3.0);
 911
 912   fitTotal->SetParameter(0,fg->GetParameter(0));
 913   fitTotal->SetParameter(1,fg->GetParameter(1));
 914   fitTotal->SetParameter(2,fg->GetParameter(2));
 915
 916   TF1* fexpo = new TF1("expo","expo",xmin,xmax);
 917
 918   hfit->Fit(fexpo,"","",3.5,5);
 919
 920   fitTotal->SetParameter(8,fexpo->GetParameter(0));
 921   fitTotal->SetParameter(9,fexpo->GetParameter(1));
 922
 923   fitTotal->SetParameter(3,h.GetMaximum()),
 924   fitTotal->SetParameter(4,cbalpha);
 925   fitTotal->SetParameter(5,cbn);
 926   fitTotal->SetParameter(6,3.15);
 927   fitTotal->SetParameter(7,0.1);
 928
 929   const char* fitOption = "SI+";
 930
 931   TFitResultPtr fitResult = hfit->Fit(fitTotal,fitOption,"",2,5);
 932
 933   r->Set("MeanJpsi",fitTotal->GetParameter(6),fitTotal->GetParError(6));
 934   r->Set("SigmaJpsi",fitTotal->GetParameter(7),fitTotal->GetParError(7));
 935
 936   double m = r->GetValue("MeanJpsi");
 937   double s = r->GetValue("SigmaJpsi");
 938   double n = 3.0;
 939
 940   TF1* fcb = new TF1("fcb",funcCB,xmin,xmax,5);
 941   fcb->SetParameters(fitTotal->GetParameter(3),
 942                      fitTotal->GetParameter(4),
 943                      fitTotal->GetParameter(5),
 944                      fitTotal->GetParameter(6),
 945                      fitTotal->GetParameter(7));
 946
 947   fcb->SetLineColor(6);
 948   fcb->SetNpx(100);
 949   TLine* l1 = new TLine(m-n*s,0,m-n*s,fitTotal->GetParameter(3));
 950   TLine* l2 = new TLine(m+n*s,0,m+n*s,fitTotal->GetParameter(3));
 951   l1->SetLineColor(6);
 952   l2->SetLineColor(6);
 953   h.GetListOfFunctions()->Add(fcb);
 954   h.GetListOfFunctions()->Add(l1);
 955   h.GetListOfFunctions()->Add(l2);
 956
 957
 958   Double_t cbParameters[5];
 959   Double_t covarianceMatrix[5][5];
 960
 961   cbParameters[0] = fitTotal->GetParameter(3);
 962   cbParameters[1] = fitTotal->GetParameter(4);
 963   cbParameters[2] = fitTotal->GetParameter(5);
 964   cbParameters[3] = fitTotal->GetParameter(6);
 965   cbParameters[4] = fitTotal->GetParameter(7);
 966
 967   for ( int iy = 0; iy < 5; ++iy )
 968   {
 969     for ( int ix = 0; ix < 5; ++ix )
 970     {
 971       covarianceMatrix[ix][iy] = (fitResult->GetCovarianceMatrix())(ix+3,iy+3);
 972     }
 973   }
 974
 975   double njpsi = fcb->Integral(m-n*s,m+n*s)/h.GetBinWidth(1);
 976
 977   double nerr = fcb->IntegralError(m-n*s,m+n*s,&cbParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(1);
 978
 979   r->Set("NofJpsi",njpsi,nerr);
 980
 981   return r;
 982 }
 983
 984 /*
 985 //_____________________________________________________________________________
 986 void AliAnalysisMuMuResult::FitJpsiPCBE(TH1& h)
 987 {
 988   std::cout << "Fit with jpsi alone (pol2 + CB + expo)" << std::endl;
 989
 990   const Double_t xmin(2.0);
 991   const Double_t xmax(5.0);
 992
 993   fitTotal = new TF1("fitTotal",funcJpsiPCBE,xmin,xmax,10);
 994   fitTotal->SetParNames("p0","p1","p2","N","alpha","n","mean","sigma","expocste","exposlope");
 995
 996   fitTotal->SetParLimits(3,0,h.GetMaximum()*2); // N
 997
 998   const Double_t cbalpha(0.98);
 999   const Double_t cbn(5.2);
1000
1001   fitTotal->FixParameter(4,cbalpha);
1002   fitTotal->FixParameter(5,cbn);
1003
1004   fitTotal->SetParLimits(6,2,4); // mean
1005   fitTotal->SetParLimits(7,0.05,0.2); // sigma
1006
1007   TF1* fpol2 = new TF1("pol2","pol2",xmin,xmax);
1008
1009   h.Fit(fpol2,"+","",2,2.8);
1010
1011   fitTotal->SetParameter(0,fpol2->GetParameter(0));
1012   fitTotal->SetParameter(1,fpol2->GetParameter(1));
1013   fitTotal->SetParameter(2,fpol2->GetParameter(2));
1014
1015   TF1* fexpo = new TF1("expo","expo",xmin,xmax);
1016
1017   h.Fit(fexpo,"+","",3.5,4.5);
1018
1019   fitTotal->SetParameter(8,fexpo->GetParameter(0));
1020   fitTotal->SetParameter(9,fexpo->GetParameter(1));
1021
1022   fitTotal->SetParameter(3,h.GetMaximum()),
1023   fitTotal->SetParameter(4,cbalpha);
1024   fitTotal->SetParameter(5,cbn);
1025   fitTotal->SetParameter(6,3.15);
1026   fitTotal->SetParameter(7,0.1);
1027
1028   h.Fit(fitTotal,"+","",2.5,5);
1029
1030   Set("MeanJpsi",fitTotal->GetParameter(6),fitTotal->GetParError(6));
1031   Set("SigmaJpsi",fitTotal->GetParameter(7),fitTotal->GetParError(7));
1032
1033   double m = GetValue("MeanJpsi");
1034   double s = GetValue("SigmaJpsi");
1035   double n = 2.0;
1036
1037   TF1* fcb = new TF1("fcb",funcCB,xmin,xmax,5);
1038   fcb->SetParameters(fitTotal->GetParameter(3),
1039                      fitTotal->GetParameter(4),
1040                      fitTotal->GetParameter(5),
1041                      fitTotal->GetParameter(6),
1042                      fitTotal->GetParameter(7));
1043
1044   fcb->SetLineColor(6);
1045   fcb->SetNpx(100);
1046   TLine* l1 = new TLine(m-n*s,0,m-n*s,fitTotal->GetParameter(3));
1047   TLine* l2 = new TLine(m+n*s,0,m+n*s,fitTotal->GetParameter(3));
1048   l1->SetLineColor(6);
1049   l2->SetLineColor(6);
1050   h.GetListOfFunctions()->Add(fcb);
1051   h.GetListOfFunctions()->Add(l1);
1052   h.GetListOfFunctions()->Add(l2);
1053
1054
1055   Set("NofJpsi",fitTotal->Integral(m-n*s,m+n*s)/h.GetBinWidth(1),fitTotal->IntegralError(m-n*s,m+n*s)/h.GetBinWidth(1));
1056
1057   //  Set("NofJpsi",fitTotal->Integral(0,10)/h.GetBinWidth(1),fitTotal->IntegralError(0,10)/h.GetBinWidth(1));
1058
1059 }
1060
1061 //_____________________________________________________________________________
1062 void AliAnalysisMuMuResult::FitJpsiCBE(TH1& h)
1063 {
1064   std::cout << "Fit with jpsi alone" << std::endl;
1065
1066   const Double_t xmin(1.5);
1067   const Double_t xmax(8.0);
1068
1069   fitTotal = new TF1("fitTotal",funcJpsiCBE,xmin,xmax,7);
1070   fitTotal->SetParNames("N","alpha","n","mean","sigma","expocste","exposlope");
1071
1072 //  fitTotal->SetParameters(h.GetMaximum(),1,5,3.0,0.07,1.5,3,1,0);
1073
1074   fitTotal->SetParameters(1,1.15,3.6,3.0,0.07,1,-1);
1075
1076   fitTotal->SetParLimits(0,0,h.GetMaximum()); // N
1077 //  fitTotal->SetParLimits(1,0.1,2); // alpha
1078   fitTotal->FixParameter(1,0.98);
1079 //  fitTotal->SetParLimits(2,0.01,5); // n
1080   fitTotal->FixParameter(2,5.2);
1081   fitTotal->SetParLimits(3,2.8,3.5); // mean
1082   fitTotal->SetParLimits(4,0.05,0.2); // sigma
1083
1084   TF1* fexpo = new TF1("expo","expo",xmin,xmax);
1085
1086   h.Fit(fexpo,"+","",2,3);
1087
1088   fitTotal->SetParameter(5,fexpo->GetParameter(0));
1089   fitTotal->SetParameter(6,fexpo->GetParameter(1));
1090
1091   h.Fit(fitTotal,"+","",2,5);
1092
1093
1094   Set("MeanJpsi",fitTotal->GetParameter(3),fitTotal->GetParError(3));
1095   Set("SigmaJpsi",fitTotal->GetParameter(4),fitTotal->GetParError(4));
1096
1097   double m = GetValue("MeanJpsi");
1098   double s = GetValue("SigmaJpsi");
1099   double n = 3.0;
1100
1101   TF1* fcb = new TF1("fcb",funcCB,xmin,xmax,5);
1102   fcb->SetParameters(fitTotal->GetParameter(0),
1103                      fitTotal->GetParameter(1),
1104                      fitTotal->GetParameter(2),
1105                      fitTotal->GetParameter(3),
1106                      fitTotal->GetParameter(4));
1107
1108   fcb->SetLineColor(6);
1109   fcb->SetNpx(1000);
1110   TLine* l1 = new TLine(m-n*s,0,m-n*s,fitTotal->GetParameter(0));
1111   TLine* l2 = new TLine(m+n*s,0,m+n*s,fitTotal->GetParameter(0));
1112   l1->SetLineColor(6);
1113   l2->SetLineColor(6);
1114   h.GetListOfFunctions()->Add(fcb);
1115   h.GetListOfFunctions()->Add(l1);
1116   h.GetListOfFunctions()->Add(l2);
1117
1118
1119   Set("NofJpsi",fitTotal->Integral(m-n*s,m+n*s)/h.GetBinWidth(1),fitTotal->IntegralError(m-n*s,m+n*s)/h.GetBinWidth(1));
1120
1121   //  Set("NofJpsi",fitTotal->Integral(0,10)/h.GetBinWidth(1),fitTotal->IntegralError(0,10)/h.GetBinWidth(1));
1122
1123 }
1124
1125 //_____________________________________________________________________________
1126 void AliAnalysisMuMuResult::FitJpsiECBE(TH1& h)
1127 {
1128   std::cout << "Fit with jpsi alone (expo + CB + expo)" << std::endl;
1129
1130   const Double_t xmin(1.5);
1131   const Double_t xmax(8.0);
1132
1133   fitTotal = new TF1("fitTotal",funcJpsiECBE,xmin,xmax,9);
1134   fitTotal->SetParNames("e0","s0","N","alpha","n","mean","sigma","e1","s1");
1135
1136   fitTotal->SetParameters(1,-1,1,1.15,3.6,3.2,0.06,-1);
1137
1138   fitTotal->SetParLimits(0,0,h.GetMaximum()*2);
1139
1140   fitTotal->FixParameter(3,0.98); // alpha
1141   fitTotal->FixParameter(4,5.2); // n
1142   fitTotal->SetParLimits(5,2.8,3.5); // mean
1143   fitTotal->SetParLimits(6,0.05,0.2); // sigma
1144
1145   TF1* fexpo1 = new TF1("expo1","expo",xmin,xmax);
1146   TF1* fexpo2 = new TF1("expo2","expo",xmin,xmax);
1147
1148   h.Fit(fexpo1,"","",1.5,3);
1149
1150   fitTotal->SetParameter(0,fexpo1->GetParameter(0));
1151   fitTotal->SetParameter(1,fexpo1->GetParameter(1));
1152
1153   h.Fit(fexpo2,"","",3.5,5.0);
1154
1155   fitTotal->SetParameter(7,fexpo2->GetParameter(0));
1156   fitTotal->SetParameter(8,fexpo2->GetParameter(1));
1157
1158   const char* fitOption = "SI+";
1159
1160   TFitResultPtr r = h.Fit(fitTotal,fitOption,"",2,5);
1161
1162   Set("MeanJpsi",fitTotal->GetParameter(5),fitTotal->GetParError(5));
1163   Set("SigmaJpsi",fitTotal->GetParameter(6),fitTotal->GetParError(6));
1164
1165   double m = GetValue("MeanJpsi");
1166   double s = GetValue("SigmaJpsi");
1167   double n = 3.0;
1168
1169   TF1* fcb = new TF1("fcb",funcCB,xmin,xmax,5);
1170   fcb->SetParameters(fitTotal->GetParameter(2),
1171                      fitTotal->GetParameter(3),
1172                      fitTotal->GetParameter(4),
1173                      fitTotal->GetParameter(5),
1174                      fitTotal->GetParameter(6));
1175
1176   fcb->SetParError(0,fitTotal->GetParError(2));
1177   fcb->SetParError(1,fitTotal->GetParError(3));
1178   fcb->SetParError(2,fitTotal->GetParError(4));
1179   fcb->SetParError(3,fitTotal->GetParError(5));
1180   fcb->SetParError(4,fitTotal->GetParError(6));
1181
1182   fcb->SetLineColor(6);
1183   fcb->SetNpx(1000);
1184   TLine* l1 = new TLine(m-n*s,0,m-n*s,fitTotal->GetParameter(2));
1185   TLine* l2 = new TLine(m+n*s,0,m+n*s,fitTotal->GetParameter(2));
1186   l1->SetLineColor(6);
1187   l2->SetLineColor(6);
1188   h.GetListOfFunctions()->Add(fcb);
1189   h.GetListOfFunctions()->Add(l1);
1190   h.GetListOfFunctions()->Add(l2);
1191
1192
1193   Double_t cbParameters[5];
1194   Double_t covarianceMatrix[5][5];
1195
1196   cbParameters[0] = fitTotal->GetParameter(2);
1197   cbParameters[1] = fitTotal->GetParameter(3);
1198   cbParameters[2] = fitTotal->GetParameter(4);
1199   cbParameters[3] = fitTotal->GetParameter(5);
1200   cbParameters[4] = fitTotal->GetParameter(6);
1201
1202   for ( int iy = 0; iy < 5; ++iy )
1203   {
1204     for ( int ix = 0; ix < 5; ++ix )
1205     {
1206       covarianceMatrix[ix][iy] = (r->GetCovarianceMatrix())(ix+2,iy+2);
1207     }
1208   }
1209
1210
1211   double njpsi = fcb->Integral(m-n*s,m+n*s)/h.GetBinWidth(1);
1212
1213   double nerr = fcb->IntegralError(m-n*s,m+n*s,&cbParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(1);
1214
1215
1216   Set("NofJpsi",njpsi,nerr);
1217 }
1218
1219  */
1220
1221 //_____________________________________________________________________________
1222 AliAnalysisMuMuResult* AliAnalysisMuMuResult::FitJpsi(TH1& h)
1223 {
1224   std::cout << "Fit with jpsi alone" << std::endl;
1225
1226   Int_t nrebin = fMinv->GetXaxis()->GetNbins() / h.GetXaxis()->GetNbins();
1227
1228   AliAnalysisMuMuResult* r = new AliAnalysisMuMuResult(h,"JPSI",nrebin);
1229
1230   TH1* hfit = r->Minv();
1231
1232   const Double_t xmin(1.5);
1233   const Double_t xmax(8.0);
1234
1235   TF1* fitTotal = new TF1("fitTotal",funcCB2,xmin,xmax,7);
1236   fitTotal->SetParNames("N","alpha","n","mean","sigma","alphaprime","nprime");
1237   fitTotal->SetParameters(h.GetMaximum(),1,5,3.0,0.07,1.5,3);
1238   fitTotal->SetParLimits(0,0,h.GetMaximum()*2); // N
1239   fitTotal->SetParLimits(1,0,10); // alpha
1240   fitTotal->SetParLimits(2,1,10); // n
1241   fitTotal->SetParLimits(3,1,4); // mean
1242   fitTotal->SetParLimits(4,0.01,1); // sigma
1243   fitTotal->SetParLimits(5,0,10); // alpha
1244   fitTotal->SetParLimits(6,1,10); // n
1245
1246   hfit->Fit(fitTotal,"+","",2,5);
1247
1248
1249   r->Set("MeanJpsi",fitTotal->GetParameter(3),fitTotal->GetParError(3));
1250   r->Set("SigmaJpsi",fitTotal->GetParameter(4),fitTotal->GetParError(4));
1251
1252   double m = r->GetValue("MeanJpsi");
1253   double s = r->GetValue("SigmaJpsi");
1254   double n = 10.0;
1255
1256   r->Set("NofJpsi",fitTotal->Integral(m-n*s,m+n*s)/h.GetBinWidth(1),fitTotal->IntegralError(m-n*s,m+n*s)/h.GetBinWidth(1));
1257
1258   return r;
1259 }
1260
1261 //_____________________________________________________________________________
1262 AliAnalysisMuMuResult* AliAnalysisMuMuResult::FitJpsi2CB2VWG(const TH1& h)
1263 {
1264   std::cout << "Fit with jpsi + psiprime VWG" << std::endl;
1265
1266   Int_t nrebin = fMinv->GetXaxis()->GetNbins() / h.GetXaxis()->GetNbins();
1267
1268   AliAnalysisMuMuResult* r = new AliAnalysisMuMuResult(h,"JPSI2CB2VWG",nrebin);
1269
1270
1271   TH1* hfit = r->Minv();
1272
1273   const Double_t xmin(2.2);
1274   const Double_t xmax(5.0);
1275
1276   TF1* fitTotal = new TF1("fitTotal",func2CB2VWG,xmin,xmax,12);
1277   fitTotal->SetParNames("kVWG","mVWG","sVWG1","sVWG2","kPsi","mPsi","sPsi","alPsi","nlPsi","auPsi","nuPsi");
1278   fitTotal->SetParName(11, "kPsi'");
1279
1280   fitTotal->SetParameter(0, 10000.);
1281   fitTotal->SetParameter(1, 1.9);
1282   fitTotal->SetParameter(2, 0.5);
1283   fitTotal->SetParLimits(2, 0., 100.);
1284   fitTotal->SetParameter(3, 0.3);
1285   fitTotal->SetParLimits(3, 0., 100.);
1286   fitTotal->SetParameter(4, 100.);
1287   fitTotal->SetParameter(5, 3.13);
1288   fitTotal->SetParLimits(5, 3.08, 3.2);
1289   fitTotal->SetParameter(6, 0.08);
1290   fitTotal->SetParLimits(6, 0.05, 0.15);
1291   fitTotal->FixParameter(7, 0.93);
1292   fitTotal->FixParameter(8, 5.59);
1293   fitTotal->FixParameter(9, 2.32);
1294   fitTotal->FixParameter(10, 3.39);
1295   fitTotal->SetParameter(11, 10.);
1296
1297   const char* fitOption = "SIER"; //+";
1298
1299   TFitResultPtr fitResult = hfit->Fit(fitTotal,fitOption,"");
1300
1301   r->Set("MeanJpsi",fitTotal->GetParameter(5),fitTotal->GetParError(5));
1302   r->Set("SigmaJpsi",fitTotal->GetParameter(6),fitTotal->GetParError(6));
1303
1304   double m = r->GetValue("MeanJpsi");
1305   double s = r->GetValue("SigmaJpsi");
1306   double n = 3.0;
1307
1308   TF1* fcb = new TF1("fcb",CrystalBallExtended,xmin,xmax,7);
1309   fcb->SetParameters(fitTotal->GetParameter(4),
1310                      fitTotal->GetParameter(5),
1311                      fitTotal->GetParameter(6),
1312                      fitTotal->GetParameter(7),
1313                      fitTotal->GetParameter(8),
1314                      fitTotal->GetParameter(9),
1315                      fitTotal->GetParameter(10));
1316
1317
1318   fcb->SetLineColor(1);
1319   fcb->SetNpx(1000);
1320   TLine* l1 = new TLine(m-n*s,0,m-n*s,fitTotal->GetParameter(4));
1321   TLine* l2 = new TLine(m+n*s,0,m+n*s,fitTotal->GetParameter(4));
1322   l1->SetLineColor(6);
1323   l2->SetLineColor(6);
1324   hfit->GetListOfFunctions()->Add(fcb);
1325   hfit->GetListOfFunctions()->Add(l1);
1326   hfit->GetListOfFunctions()->Add(l2);
1327
1328   Double_t cbParameters[7];
1329   Double_t covarianceMatrix[7][7];
1330
1331   for ( int ix = 0; ix < 7; ++ix )
1332   {
1333     cbParameters[ix] = fitTotal->GetParameter(ix+4);
1334   }
1335
1336   for ( int iy = 0; iy < 5; ++iy )
1337   {
1338     for ( int ix = 0; ix < 5; ++ix )
1339     {
1340       covarianceMatrix[ix][iy] = (fitResult->GetCovarianceMatrix())(ix+4,iy+4);
1341     }
1342   }
1343
1344   double njpsi = fcb->Integral(m-n*s,m+n*s)/h.GetBinWidth(1);
1345   double nerr = fcb->IntegralError(m-n*s,m+n*s,&cbParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(1);
1346
1347   r->Set("NofJpsi",njpsi,nerr);
1348
1349   return r;
1350 }
1351
1352 /*
1353 //_____________________________________________________________________________
1354 void AliAnalysisMuMuResult::FitUpsilon(TH1& h)
1355 {
1356   std::cout << "Fit with upsilon alone" << std::endl;
1357
1358   const Double_t xmin(6.0);
1359   const Double_t xmax(12.0);
1360
1361   fitTotal = new TF1("fitTotal",funcCB2,xmin,xmax,7);
1362   fitTotal->SetParNames("N","alpha","n","mean","sigma","alphaprime","nprime");
1363   fitTotal->SetParameters(h.GetMaximum(),1,5,9.46,0.2,1.5,3);
1364   fitTotal->SetParLimits(0,0,h.GetMaximum()*2); // N
1365   fitTotal->SetParLimits(1,0,10); // alpha
1366   fitTotal->SetParLimits(2,1,10); // n
1367   fitTotal->SetParLimits(3,8,12); // mean
1368   fitTotal->SetParLimits(4,0.01,1); // sigma
1369   fitTotal->SetParLimits(5,0,10); // alpha
1370   fitTotal->SetParLimits(6,1,10); // n
1371
1372   h.Fit(fitTotal,"+","",6,12);
1373
1374
1375   Set("MeanUpsilon",fitTotal->GetParameter(3),fitTotal->GetParError(3));
1376   Set("SigmaUpsilon",fitTotal->GetParameter(4),fitTotal->GetParError(4));
1377
1378   double m = GetValue("MeanUpsilon");
1379   double s = GetValue("SigmaUpsilon");
1380   double n = 3.0;
1381
1382   Set("NofUpsilon",fitTotal->Integral(m-n*s,m+n*s)/h.GetBinWidth(1),fitTotal->IntegralError(m-n*s,m+n*s)/h.GetBinWidth(1));
1383 }
1384 */
1385
1386 //_____________________________________________________________________________
1387 Double_t AliAnalysisMuMuResult::ErrorAB(Double_t a, Double_t aerr, Double_t b, Double_t berr)
1388 {
1389   Double_t e(0.0);
1390
1391   if ( TMath::Abs(a) > 1E-12 )
1392   {
1393     e += (aerr*aerr)/(a*a);
1394   }
1395
1396   if ( TMath::Abs(b) > 1E-12 )
1397   {
1398     e += (berr*berr)/(b*b);
1399   }
1400
1401   return TMath::Sqrt(e);
1402 }
1403
1404 //_____________________________________________________________________________
1405 Double_t AliAnalysisMuMuResult::ErrorABC(Double_t a, Double_t aerr, Double_t b, Double_t berr, Double_t c, Double_t cerror)
1406 {
1407   Double_t e(0.0);
1408
1409   if ( TMath::Abs(a) > 1E-12 )
1410   {
1411     e += (aerr*aerr)/(a*a);
1412   }
1413
1414   if ( TMath::Abs(b) > 1E-12 )
1415   {
1416     e += (berr*berr)/(b*b);
1417   }
1418
1419   if ( TMath::Abs(b) > 1E-12 )
1420   {
1421     e += (cerror*cerror)/(c*c);
1422   }
1423
1424   return TMath::Sqrt(e);
1425 }
1426
1427
1428 //_____________________________________________________________________________
1429 Bool_t AliAnalysisMuMuResult::AddFit(const char* fitType, Int_t nrebin)
1430 {
1431   //  new TCanvas;
1432
1433   if ( !fMinv ) return kFALSE;
1434
1435   TString sFitType(fitType);
1436   sFitType.ToUpper();
1437
1438   if ( fMinv->GetEntries()<100 && !sFitType.Contains("COUNT")) return kFALSE;
1439
1440   TH1* hminv = static_cast<TH1*>(fMinv->Clone());
1441
1442   hminv->Rebin(nrebin);
1443   hminv->SetDirectory(0);
1444
1445   AliAnalysisMuMuResult* r(0x0);
1446
1447 //  if ( sFitType == "PSI12" )
1448 //  {
1449 //    r = FitJpsiPsiPrimeCustom(*hminv);
1450 //  }
1451 //
1452   if ( sFitType=="PSI1")
1453   {//
1454     r = FitJpsi(*hminv);
1455 //    r = FitJpsiGCBE(*hminv);
1456   }
1457
1458   if ( sFitType == "PSILOW")
1459   {
1460     r = FitJpsi2CB2VWG(*hminv);
1461   }
1462
1463   if ( sFitType == "COUNTPSI" )
1464   {
1465     r = new AliAnalysisMuMuResult(*hminv,"COUNTJPSI",1);
1466     Double_t n = CountParticle(*hminv,"Jpsi");
1467     r->Set("NofJpsi",n,TMath::Sqrt(n));
1468   }
1469
1470 //  if (r)
1471 //  {
1472 //    r->SetNofInputParticles("Jpsi",GetValue("NofInputJpsi"));
1473 //  }
1474
1475   if (!fSubResults)
1476   {
1477     fSubResults = new TObjArray;
1478     fSubResults->SetOwner(kTRUE);
1479   }
1480
1481   if ( r )
1482   {
1483     r->Print();
1484     r->SetBin(Bin());
1485     r->SetNofTriggers(NofTriggers());
1486     r->SetNofRuns(NofRuns());
1487     fSubResults->Add(r);
1488   }
1489
1490   delete hminv;
1491
1492   return (r!=0x0);
1493 }
1494
1495 //_____________________________________________________________________________
1496 Double_t AliAnalysisMuMuResult::GetErrorStat(const char* name, const char* subResultName) const
1497 {
1498   // compute the mean value from all subresults
1499
1500   if ( strlen(subResultName) > 0 )
1501   {
1502     if ( !fSubResults)
1503     {
1504       AliError(Form("No subresult from which I could get the %s one...",subResultName));
1505       return TMath::Limits<Double_t>::Max();
1506     }
1507     AliAnalysisMuMuResult* sub = static_cast<AliAnalysisMuMuResult*>(fSubResults->FindObject(subResultName));
1508     if (!sub)
1509     {
1510       AliError(Form("Could not get subresult named %s",subResultName));
1511       return TMath::Limits<Double_t>::Max();
1512     }
1513     return sub->GetErrorStat(name);
1514   }
1515
1516   if ( fMap )
1517   {
1518     TObjArray* p = static_cast<TObjArray*>(fMap->GetValue(name));
1519     if (p)
1520     {
1521       TParameter<double>* val = static_cast<TParameter<double>*>(p->At(kErrorStat));
1522       return val->GetVal();
1523     }
1524   }
1525
1526   TIter next(fSubResults);
1527   AliAnalysisMuMuResult* r;
1528   Int_t n(0);
1529   Double_t mean(0);
1530
1531   while ( ( r = static_cast<AliAnalysisMuMuResult*>(next()) ) )
1532   {
1533     if ( r->HasValue(name) )
1534     {
1535       mean += r->GetErrorStat(name);
1536       ++n;
1537     }
1538   }
1539   return ( n ? mean/n : 0.0 );
1540 }
1541
1542 //_____________________________________________________________________________
1543 Double_t AliAnalysisMuMuResult::GetValue(const char* name, const char* subResultName) const
1544 {
1545   // get a value (either directly or by computing the mean of the subresults)
1546
1547   if ( strlen(subResultName) > 0 )
1548   {
1549     if ( !fSubResults)
1550     {
1551       AliError(Form("No subresult from which I could get the %s one...",subResultName));
1552       return TMath::Limits<Double_t>::Max();
1553     }
1554     AliAnalysisMuMuResult* sub = static_cast<AliAnalysisMuMuResult*>(fSubResults->FindObject(subResultName));
1555     if (!sub)
1556     {
1557       AliError(Form("Could not get subresult named %s",subResultName));
1558       return TMath::Limits<Double_t>::Max();
1559     }
1560     return sub->GetValue(name);
1561   }
1562
1563   if (fMap)
1564   {
1565     TObjArray* p = static_cast<TObjArray*>(fMap->GetValue(name));
1566     if (p)
1567     {
1568       TParameter<double>* val = static_cast<TParameter<double>*>(p->At(kValue));
1569       return val->GetVal();
1570     }
1571   }
1572
1573   // compute the mean value from all subresults
1574   TIter next(fSubResults);
1575   AliAnalysisMuMuResult* r;
1576   Int_t n(0);
1577   Double_t mean(0);
1578
1579   while ( ( r = static_cast<AliAnalysisMuMuResult*>(next()) ) )
1580   {
1581     if ( r->HasValue(name) )
1582     {
1583       mean += r->GetValue(name);
1584       ++n;
1585     }
1586   }
1587   return ( n ? mean/n : 0.0 );
1588 }
1589
1590 //_____________________________________________________________________________
1591 Bool_t AliAnalysisMuMuResult::HasValue(const char* name, const char* subResultName) const
1592 {
1593   if ( strlen(subResultName) > 0 )
1594   {
1595     if ( !fSubResults)
1596     {
1597       AliError(Form("No subresult from which I could get the %s one...",subResultName));
1598       return kFALSE;
1599     }
1600     AliAnalysisMuMuResult* sub = static_cast<AliAnalysisMuMuResult*>(fSubResults->FindObject(subResultName));
1601     if (!sub)
1602     {
1603       AliError(Form("Could not get subresult named %s",subResultName));
1604       return kFALSE;
1605     }
1606     return sub->HasValue(name);
1607   }
1608
1609   if ( fMap && ( fMap->GetValue(name) != 0x0 ) )
1610   {
1611     return kTRUE;
1612   }
1613
1614   TIter next(fSubResults);
1615   AliAnalysisMuMuResult* r;
1616
1617   while ( ( r = static_cast<AliAnalysisMuMuResult*>(next()) ) )
1618   {
1619     if ( r->HasValue(name) ) return kTRUE;
1620   }
1621
1622   return kFALSE;
1623 }
1624
1625 //_____________________________________________________________________________
1626 THashList* AliAnalysisMuMuResult::Keys() const
1627 {
1628   /// Return the complete list of keys we're using
1629   if (!fKeys)
1630   {
1631     fKeys = new THashList;
1632     fKeys->SetOwner(kTRUE);
1633     TIter next(fMap);
1634     TObjString* key;
1635
1636     while ( ( key = static_cast<TObjString*>(next()) ) )
1637     {
1638       if ( !fKeys->FindObject(key->String()) )
1639       {
1640         fKeys->Add(new TObjString(key->String()));
1641       }
1642     }
1643
1644     AliAnalysisMuMuResult* r;
1645     TIter nextResult(fSubResults);
1646
1647     while ( ( r = static_cast<AliAnalysisMuMuResult*>(nextResult())) )
1648     {
1649       TIter nextHL(r->Keys());
1650       TObjString* s;
1651
1652       while ( ( s = static_cast<TObjString*>(nextHL())) )
1653       {
1654         if ( !fKeys->FindObject(s->String()) )
1655         {
1656           fKeys->Add(new TObjString(s->String()));
1657         }
1658       }
1659     }
1660
1661   }
1662   return fKeys;
1663 }
1664
1665 //_____________________________________________________________________________
1666 Long64_t AliAnalysisMuMuResult::Merge(TCollection* list)
1667 {
1668   /// Merge method
1669   ///
1670   /// Merge a list of AliAnalysisMuMuResult objects with this
1671   /// Returns the number of merged objects (including this).
1672   ///
1673   /// Note that the merging is to be understood here as an average operation
1674
1675   AliInfo(Form("this=%p",this));
1676   if (!list) return 0;
1677
1678   if (list->IsEmpty()) return 1;
1679
1680   TIter next(list);
1681   TObject* currObj;
1682   TList hminvList;
1683   hminvList.SetOwner(kTRUE);
1684
1685   Double_t thisWeight = Weight();
1686   Double_t sumOfWeights = thisWeight;
1687
1688   Double_t nofTriggers = fNofTriggers*thisWeight;
1689   Double_t nofRuns = fNofRuns*thisWeight;
1690
1691   fNofRuns *= thisWeight;
1692
1693   while ( ( currObj = next() ) )
1694   {
1695     AliAnalysisMuMuResult* result = dynamic_cast<AliAnalysisMuMuResult*>(currObj);
1696     if (!result)
1697     {
1698       AliFatal(Form("object named \"%s\" is a %s instead of an AliAnalysisMuMuResult!", currObj->GetName(), currObj->ClassName()));
1699       continue;
1700     }
1701
1702     Double_t w = result->Weight();
1703
1704     nofRuns += result->NofRuns()*w;
1705     nofTriggers += result->NofTriggers()*w;
1706     fWeight += result->fWeight;
1707     sumOfWeights += w;
1708   }
1709
1710   thisWeight/= sumOfWeights;
1711   fNofRuns = nofRuns/sumOfWeights;
1712   fNofTriggers = nofTriggers/sumOfWeights;
1713   fWeight /= sumOfWeights;
1714
1715   AliInfo(Form("thisWeight=%e sumOfWeight=%8.2f noftriggers=%e weight=%e",thisWeight,sumOfWeights,1.0*fNofTriggers,fWeight));
1716
1717   TIter nextKey(fMap);
1718   TObjString* key;
1719
1720   while ( ( key = static_cast<TObjString*>(nextKey())) )
1721   {
1722     AliInfo(key->String().Data());
1723
1724     Double_t value = GetValue(key->String())*thisWeight;
1725     Double_t estat = GetErrorStat(key->String())*GetErrorStat(key->String())*thisWeight*thisWeight;
1726
1727     Double_t test(thisWeight);
1728
1729     next.Reset();
1730
1731     while ( ( currObj = next() ) )
1732     {
1733       AliAnalysisMuMuResult* result = dynamic_cast<AliAnalysisMuMuResult*>(currObj);
1734
1735       if (!result->HasValue(key->String()))
1736       {
1737         AliError(Form("Did not find key %s in of the result to merge",key->String().Data()));
1738         continue;
1739       }
1740
1741       // can only merge under the condition we have the same bin
1742
1743       if ( fBin != result->Bin() )
1744       {
1745         AliError("Cannot merge results of different bin");
1746         continue;
1747       }
1748
1749       Double_t w = result->Weight()/sumOfWeights;
1750
1751       Double_t w2 = w*w;
1752
1753       test += w;
1754
1755       value += result->GetValue(key->String())*w;
1756       estat += result->GetErrorStat(key->String())*result->GetErrorStat(key->String())*w2;
1757
1758     }
1759
1760     Set(key->String(),
1761         value,
1762         TMath::Sqrt(estat));
1763
1764     AliInfo(Form("test=%e",test));
1765   }
1766
1767   if ( fMinv )
1768   {
1769     fMinv->Scale(thisWeight);
1770   }
1771
1772   next.Reset();
1773
1774   while ( ( currObj = next() ) )
1775   {
1776     AliAnalysisMuMuResult* result = dynamic_cast<AliAnalysisMuMuResult*>(currObj);
1777
1778     if ( result->Minv() )
1779     {
1780       TH1* h = static_cast<TH1*>(result->Minv()->Clone());
1781       AliInfo(Form("Nbins %d xmin %e xmax %e",h->GetXaxis()->GetNbins(),h->GetXaxis()->GetXmin(),
1782                    h->GetXaxis()->GetXmax()));
1783       h->Scale(result->Weight());
1784       hminvList.Add(h);
1785     }
1786   }
1787
1788   if ( fMinv )
1789   {
1790     fMinv->Merge(&hminvList);
1791     fMinv->Scale(1.0/sumOfWeights);
1792   }
1793
1794   TIter nextSubresult(fSubResults);
1795   AliAnalysisMuMuResult* r;
1796
1797   while ( ( r = static_cast<AliAnalysisMuMuResult*>(nextSubresult())) )
1798   {
1799     TList sublist;
1800
1801     next.Reset();
1802
1803     while ( ( currObj = next() ) )
1804     {
1805       sublist.Add(currObj);
1806     }
1807
1808     r->Merge(&sublist);
1809   }
1810
1811   return list->GetEntries()+1;
1812 }
1813
1814 //_____________________________________________________________________________
1815 Int_t AliAnalysisMuMuResult::NofRuns() const
1816 {
1817   if ( !Mother() ) return fNofRuns;
1818   else return Mother()->NofRuns();
1819 }
1820
1821 //_____________________________________________________________________________
1822 Int_t AliAnalysisMuMuResult::NofTriggers() const
1823 {
1824   if ( !Mother() ) return fNofTriggers;
1825   else return Mother()->NofTriggers();
1826 }
1827
1828 //_____________________________________________________________________________
1829 void AliAnalysisMuMuResult::Print(Option_t* opt) const
1830 {
1831   TString sopt(opt);
1832   sopt.ToUpper();
1833
1834   for ( Int_t i = 0; i < 9; ++i )
1835   {
1836     sopt.ReplaceAll(Form("%d",i),"");
1837   }
1838
1839   TString pot(sopt);
1840   pot.ReplaceAll("ALL","");
1841   pot.ReplaceAll("FULL","");
1842
1843   std::cout << pot.Data() << Form("%50s- NRUNS %d - NTRIGGER %10d - %s%s",
1844                GetName(),
1845                     NofRuns(),
1846                     NofTriggers(),
1847                                   fWeight > 0.0 ? Form(" WEIGHT %e -",fWeight) : "",
1848                            fBin.AsString().Data());
1849
1850     if ( fSubResults && fSubResults->GetEntries()>1 )
1851     {
1852       std::cout << " (" << fSubResults->GetEntries()-1 << " subresults)";
1853     }
1854
1855   std::cout << std::endl;
1856
1857   if ( sopt.Contains("DUMP") )
1858   {
1859     TIter next(fMap);
1860     TObjString* str;
1861     while ( ( str = static_cast<TObjString*>(next()) ) )
1862     {
1863       TObjArray* a = static_cast<TObjArray*>(fMap->GetValue(str->String()));
1864
1865       std::cout << Form("%s %e %e",
1866                         str->String().Data(),
1867                         static_cast<TParameter<Double_t>*> (a->At(kValue))->GetVal(),
1868                         static_cast<TParameter<Double_t>*> (a->At(kErrorStat))->GetVal()) << std::endl;
1869     }
1870   }
1871   else
1872   {
1873
1874     PrintParticle("Jpsi",pot.Data());
1875     PrintParticle("PsiPrime",pot.Data());
1876     PrintParticle("Upsilon",pot.Data());
1877   }
1878
1879   if ( HasValue("MBR"))
1880   {
1881     std::cout << Form("\t\tMBR %e +- %e",GetValue("MBR"),GetErrorStat("MBR")) << std::endl;
1882   }
1883
1884   if ( fSubResults && fSubResults->GetEntries() > 1 && ( sopt.Contains("ALL") || sopt.Contains("FULL") ) )
1885   {
1886     std::cout << pot.Data() << "\t===== sub results ===== " << std::endl;
1887
1888     sopt += "\t\t";
1889
1890     TIter next(fSubResults);
1891     AliAnalysisMuMuResult* r;
1892
1893     while ( ( r = static_cast<AliAnalysisMuMuResult*>(next()) ) )
1894     {
1895       r->Print(sopt.Data());
1896     }
1897   }
1898 }
1899
1900 //_____________________________________________________________________________
1901 void AliAnalysisMuMuResult::PrintParticle(const char* particle, const char* opt) const
1902 {
1903   Double_t npart = GetValue(Form("Nof%s",particle));
1904   if (npart<=0) return;
1905
1906
1907   std::cout << opt << Form("\t%s",particle) << std::endl;
1908
1909   //  Double_t npartError = GetErrorStat(Form("Nof%s",particle));
1910 //  std::cout << opt << Form("\t\t%20s %9.2f +- %5.2f","Count",npart,npartError) << std::endl;
1911
1912   TIter next(Keys());
1913   TObjString* key;
1914
1915   while ( ( key = static_cast<TObjString*>(next()) ) )
1916   {
1917     if ( !key->String().Contains(particle) ) continue;
1918
1919     PrintValue(key->String(),opt,GetValue(key->String()),GetErrorStat(key->String()));
1920   }
1921 }
1922
1923 //_____________________________________________________________________________
1924 void AliAnalysisMuMuResult::PrintValue(const char* key, const char* opt, Double_t value, Double_t errorStat) const
1925 {
1926   // print one value and its associated error
1927
1928   if ( TString(key).Contains("AccEff") )
1929   {
1930     std::cout << opt << Form("\t\t%20s %9.2f +- %5.2f %%",key,value*100,errorStat*100) << std::endl;
1931   }
1932   else if ( TString(key).BeginsWith("Sigma") ||  TString(key).BeginsWith("Mean") )
1933   {
1934     std::cout << opt << Form("\t\t%20s %9.2f +- %5.2f MeV/c^2",key,value*1E3,1E3*errorStat) << std::endl;
1935   }
1936   else if ( TString(key).BeginsWith("Nof") )
1937   {
1938     std::cout << opt << Form("\t\t%20s %9.2f +- %5.2f",key,value,errorStat) << std::endl;
1939   }
1940   else if ( value > 1E-3 && value < 1E3 )
1941   {
1942     std::cout << opt << Form("\t\t%20s %9.2f +- %5.2f",key,value,errorStat) << std::endl;
1943   }
1944   else
1945   {
1946     std::cout << opt << Form("\t\t%20s %9.2e +- %9.2e",key,value,errorStat) << std::endl;
1947   }
1948 }
1949
1950 //_____________________________________________________________________________
1951 void AliAnalysisMuMuResult::Set(const char* name, Double_t value, Double_t errorStat)
1952 {
1953   if ( !fMap )
1954   {
1955     fMap = new TMap;
1956     fMap->SetOwnerKeyValue(kTRUE,kTRUE);
1957   }
1958
1959   TObjArray* p = static_cast<TObjArray*>(fMap->GetValue(name));
1960   if (!p)
1961   {
1962     p = new TObjArray(4);
1963
1964     p->SetOwner(kTRUE);
1965
1966     p->AddAt(new TParameter<Double_t>(name,value),kValue);
1967     p->AddAt(new TParameter<Double_t>(name,errorStat),kErrorStat);
1968
1969     fMap->Add(new TObjString(name),p);
1970   }
1971   else
1972   {
1973     static_cast<TParameter<double>*>(p->At(kValue))->SetVal(value);
1974     static_cast<TParameter<double>*>(p->At(kErrorStat))->SetVal(errorStat);
1975   }
1976
1977   if ( TString(name)=="NofJpsi" )
1978   {
1979     if ( NofTriggers() > 0 )
1980     {
1981       Double_t rate = value/NofTriggers();
1982       Double_t rateError = rate*ErrorAB(value,errorStat,NofTriggers(),TMath::Sqrt(NofTriggers()));
1983       Set("RateJpsi",rate,rateError);
1984     }
1985   }
1986 }
1987
1988 //_____________________________________________________________________________
1989 void AliAnalysisMuMuResult::SetBin(const AliAnalysisMuMuBinning::Range& bin)
1990 {
1991   if (!Mother()) fBin = bin;
1992   else Mother()->SetBin(bin);
1993 }
1994
1995 //_____________________________________________________________________________
1996 void AliAnalysisMuMuResult::SetNofInputParticles(const TH1& hminv)
1997 {
1998   static const char* particleNames[] = { "Jpsi" , "PsiPrime", "Upsilon","UpsilonPrime" };
1999
2000   for ( Int_t i = 0; i < 4; ++i )
2001   {
2002     Double_t n = CountParticle(hminv,particleNames[i]);
2003
2004     AliDebug(1,Form("i=%d particle %s n %e",i,particleNames[i],n));
2005
2006     if ( n > 0 )
2007     {
2008       SetNofInputParticles(particleNames[i],n);
2009     }
2010   }
2011 }
2012
2013 //_____________________________________________________________________________
2014 void AliAnalysisMuMuResult::SetNofInputParticles(const char* particle, int n)
2015 {
2016   /// Set the number of input particles (so it is a MC result)
2017   /// and (re)compute the AccxEff values
2018
2019   Set(Form("NofInput%s",particle),n,TMath::Sqrt(n));
2020
2021   if (n<=0)
2022   {
2023     Set(Form("AccEff%s",particle),0,0);
2024     return;
2025   }
2026
2027   Double_t npart = GetValue(Form("Nof%s",particle));
2028   Double_t npartErr  = GetErrorStat(Form("Nof%s",particle));
2029   Double_t ninput = GetValue(Form("NofInput%s",particle));
2030   Double_t ninputErr = GetErrorStat(Form("NofInput%s",particle));
2031
2032   Set(Form("AccEff%s",particle),
2033       npart/ninput,
2034       (npart/ninput)*ErrorAB(npart,npartErr,ninput,ninputErr));
2035
2036   TIter next(fSubResults);
2037   AliAnalysisMuMuResult* r;
2038
2039   while ( ( r = static_cast<AliAnalysisMuMuResult*>(next())) )
2040   {
2041     r->Set(Form("NofInput%s",particle),n,TMath::Sqrt(n));
2042
2043     npart = r->GetValue(Form("Nof%s",particle));
2044     npartErr = r->GetErrorStat(Form("Nof%s",particle));
2045
2046     r->Set(Form("AccEff%s",particle),
2047            npart/ninput,
2048            (npart/ninput)*ErrorAB(npart,npartErr,ninput,ninputErr));
2049
2050   }
2051 }
2052
2053 //_____________________________________________________________________________
2054 void AliAnalysisMuMuResult::SetNofRuns(Int_t n)
2055 {
2056   if ( !Mother() ) fNofRuns=n;
2057   else Mother()->SetNofRuns(n);
2058 }
2059
2060 //_____________________________________________________________________________
2061 void AliAnalysisMuMuResult::SetNofTriggers(Int_t n)
2062 {
2063   if ( !Mother() ) fNofTriggers=n;
2064   else Mother()->SetNofTriggers(n);
2065 }
2066
2067 //_____________________________________________________________________________
2068 void AliAnalysisMuMuResult::SetMinv(const TH1& hminv)
2069 {
2070     /// Set the inv. mass spectrum to be fitted.
2071   static UInt_t n(0);
2072
2073   delete fMinv;
2074   fMinv = static_cast<TH1*>(hminv.Clone(Form("Minv%u",n++)));
2075   fMinv->SetDirectory(0);
2076 }
2077