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   /// a simple map of masses...
  42   static std::map<std::string,Double_t> massMap;
  43   // not super elegant, but that way we do not depend on TDatabasePDG and thus
  44   // can decide on our particle naming
  45   if (massMap.empty())
  46   {
  47     massMap["Jpsi"] = 3.096916e+00;
  48     massMap["PsiPrime"] = 3.68609e+00;
  49     massMap["Upsilon"] = 9.46030e+00;
  50     massMap["UpsilonPrime"] = 1.00233e+01;
  51   }
  52   return massMap;
  53 }
  54
  55
  56 Double_t funcCB(Double_t* xx, Double_t* par)
  57 {
  58   /// Crystal ball
  59
  60   Double_t norm = par[0];
  61   Double_t alpha = par[1];
  62   Double_t n = par[2];
  63   Double_t mean = par[3];
  64   Double_t sigma = par[4];
  65
  66   Double_t x = xx[0];
  67
  68   Double_t a = TMath::Power(n/TMath::Abs(alpha),n)*TMath::Exp(-0.5*alpha*alpha);
  69   Double_t b = n/TMath::Abs(alpha) - TMath::Abs(alpha);
  70
  71   Double_t y = ( TMath::Abs(sigma) > 1E-12 ? (x-mean)/sigma : 0 );
  72
  73   if ( y > alpha*-1.0 )
  74   {
  75     return norm*TMath::Exp(-0.5*y*y);
  76   }
  77   else
  78   {
  79     return norm*a*TMath::Power(b-y,-n);
  80   }
  81 }
  82
  83 Double_t funcJpsiGCBE(Double_t* xx, Double_t* par)
  84 {
  85   /// crystal ball + expo + gaussian
  86   Double_t x = xx[0];
  87
  88   Double_t g = par[0]*TMath::Gaus(x,par[1],par[2]);
  89
  90   Double_t jpsi = funcCB(xx,par+3);
  91
  92   Double_t expo = par[8]*TMath::Exp(par[9]*x);
  93
  94   return g+expo+jpsi;
  95 }
  96
  97 Double_t funcJpsiPsiPrimeCustom(Double_t* xx, Double_t* par)
  98 {
  99   // custom fit for jpsi + psi prime
 100
 101   Double_t norm = par[0];
 102   Double_t alpha = par[1];
 103   Double_t n = par[2];
 104   Double_t mean = par[3];
 105   Double_t sigma = par[4];
 106   Double_t alphaprime = par[5];
 107   Double_t nprime = par[6];
 108
 109   Double_t x = xx[0];
 110
 111   Double_t a = TMath::Power(n/TMath::Abs(alpha),n)*TMath::Exp(-0.5*alpha*alpha);
 112   Double_t b = n/TMath::Abs(alpha) - TMath::Abs(alpha);
 113   Double_t c = TMath::Power(nprime/TMath::Abs(alphaprime),nprime)*TMath::Exp(-0.5*alphaprime*alphaprime);
 114   Double_t d = nprime/TMath::Abs(alphaprime) - TMath::Abs(alphaprime);
 115
 116   Double_t y = ( TMath::Abs(sigma) > 1E-12 ? (x-mean)/sigma : 0 );
 117
 118   Double_t cb(0);
 119
 120   if ( y > alphaprime )
 121   {
 122     cb = norm*c*TMath::Power(d+y,-nprime);
 123   }
 124   else if ( y > alpha*-1.0 )
 125   {
 126     cb = norm*TMath::Exp(-0.5*y*y);
 127   }
 128   else
 129   {
 130     cb = norm*a*TMath::Power(b-y,-n);
 131   }
 132
 133   if ( x < mean )
 134   {
 135     return cb + par[7] + par[8]*x; // gaus + pol1
 136   }
 137   else
 138   {
 139     Double_t yprime = (x-par[10])/par[11];
 140     return cb + par[9]*TMath::Exp(-0.5*yprime*yprime)+par[12]*TMath::Exp(-par[13]*x);
 141     // gaus (j/psi) + gaus (psi') + expo
 142   }
 143 }
 144
 145
 146 Double_t funcCB2(Double_t* xx, Double_t* par)
 147 {
 148   /// CB2 = extended crystal ball
 149
 150   Double_t norm = par[0];
 151   Double_t alpha = par[1];
 152   Double_t n = par[2];
 153   Double_t mean = par[3];
 154   Double_t sigma = par[4];
 155   Double_t alphaprime = par[5];
 156   Double_t nprime = par[6];
 157
 158   Double_t x = xx[0];
 159
 160   Double_t a = TMath::Power(n/TMath::Abs(alpha),n)*TMath::Exp(-0.5*alpha*alpha);
 161   Double_t b = n/TMath::Abs(alpha) - TMath::Abs(alpha);
 162   Double_t c = TMath::Power(nprime/TMath::Abs(alphaprime),nprime)*TMath::Exp(-0.5*alphaprime*alphaprime);
 163   Double_t d = nprime/TMath::Abs(alphaprime) - TMath::Abs(alphaprime);
 164
 165   Double_t y = ( TMath::Abs(sigma) > 1E-12 ? (x-mean)/sigma : 0 );
 166
 167   if ( y > alphaprime )
 168   {
 169     return norm*c*TMath::Power(d+y,-nprime);
 170   }
 171   else if ( y > alpha*-1.0 )
 172   {
 173     return norm*TMath::Exp(-0.5*y*y);
 174   }
 175   else
 176   {
 177     return norm*a*TMath::Power(b-y,-n);
 178   }
 179 }
 180
 181 Double_t funcJpsiPsiPrime(Double_t* xx, Double_t* par)
 182 {
 183   /// CB + CB2
 184
 185   Double_t jpsi = funcCB(xx,par);
 186   Double_t psiprime = funcCB2(xx,par+5);
 187
 188   int n = 10;
 189   Double_t x = xx[0];
 190
 191   Double_t e1 = par[n]*TMath::Exp(par[n+1]*x);
 192   Double_t e2 = par[n+2]*TMath::Exp(par[n+3]*x);
 193
 194   Double_t e = e1;
 195
 196   if ( x > par[3] ) e=e2;
 197
 198   return jpsi+psiprime+e;
 199 }
 200
 201 Double_t funcJpsiCBE(Double_t* xx, Double_t* par)
 202 {
 203   // CB + expo
 204
 205   Double_t jpsi = funcCB(xx,par);
 206
 207   Double_t x = xx[0];
 208
 209   Double_t e1 = par[5]*TMath::Exp(par[6]*x);
 210
 211   return jpsi+e1;
 212 }
 213
 214
 215 Double_t funcJpsiPCBE(Double_t* xx, Double_t* par)
 216 {
 217   // CB + expo + pol2
 218
 219   Double_t x = xx[0];
 220
 221   Double_t pol2 = par[0] + par[1]*x + par[2]*x*x;
 222
 223   Double_t jpsi = funcCB(xx,par+3);
 224
 225   Double_t expo = par[8]*TMath::Exp(par[9]*x);
 226
 227   return pol2+jpsi+expo;
 228 }
 229
 230 Double_t funcJpsiECBE(Double_t* xx, Double_t* par)
 231 {
 232   // CB + expo
 233
 234   Double_t jpsi = funcCB(xx,par+2);
 235
 236   Double_t x = xx[0];
 237
 238   Double_t e1 = par[0]*TMath::Exp(par[1]*x);
 239
 240   Double_t e2 = par[7]*TMath::Exp(par[8]*x);
 241
 242   return e1+e2+jpsi;
 243 }
 244
 245 Double_t funcJpsiNA48(Double_t*x, Double_t* par)
 246 {
 247   /// Fit function from e.q. 4.8 of Ruben's PhD.
 248   Double_t c1 = par[0];
 249   Double_t c2 = par[1];
 250   Double_t d1 = par[2];
 251   Double_t d2 = par[3];
 252   Double_t g1 = par[4];
 253   Double_t g2 = par[5];
 254   Double_t m0 = par[6];
 255   Double_t sigma1 = par[7];
 256   Double_t sigma2 = par[8];
 257   Double_t b1 = par[9];
 258   Double_t b2 = par[10];
 259   Double_t norm = par[11];
 260
 261   Double_t m = x[0];
 262
 263   Double_t rv(0);
 264
 265   if ( m <= c1*m0 )
 266   {
 267     Double_t e = d1-g1*TMath::Sqrt(c1*m0-m);
 268     rv = TMath::Power(b1*(c1*m0-m),e);
 269     rv += sigma1;
 270   }
 271   else if( m >= c1*m0 && m <= m0 )
 272   {
 273     rv = sigma1;
 274   }
 275   else if ( m >= m0 && m < c2*m0 )
 276   {
 277     rv = sigma2;
 278   }
 279   else if( m >= c2*m0 )
 280   {
 281     Double_t e = d2-g2*TMath::Sqrt(m-c2*m0);
 282     rv = TMath::Power(b2*(m-c2*m0),e);
 283     rv += sigma2;
 284   }
 285
 286   return norm*TMath::Exp(-(m-m0)*(m-m0)/(2.0*rv*rv));
 287 }
 288
 289 const char* NormalizeName(const char* name, const char* suffix)
 290 {
 291   /// Remove - and / from the name, and adds _suffix
 292
 293   TString str(Form("%s_%s",name,suffix));
 294
 295   str.ReplaceAll("-","_");
 296   str.ReplaceAll("/","%");
 297
 298   return str.Data();
 299 }
 300
 301 //------------------------------------------------------------------------------
 302 Double_t BackgroundVWG(Double_t *x, Double_t *par)
 303 {
 304   // gaussian variable width
 305   Double_t sigma = par[2]+par[3]*((x[0]-par[1])/par[1]);
 306   return par[0]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/(2.*sigma*sigma));
 307
 308 }
 309
 310 //------------------------------------------------------------------------------
 311 Double_t CrystalBallExtended(Double_t *x,Double_t *par)
 312 {
 313   //par[0] = Normalization
 314   //par[1] = mean
 315   //par[2] = sigma
 316   //par[3] = alpha
 317   //par[4] = n
 318   //par[5] = alpha'
 319   //par[6] = n'
 320
 321   Double_t t = (x[0]-par[1])/par[2];
 322   if (par[3] < 0) t = -t;
 323
 324   Double_t absAlpha = fabs((Double_t)par[3]);
 325   Double_t absAlpha2 = fabs((Double_t)par[5]);
 326
 327   if (t >= -absAlpha && t < absAlpha2) // gaussian core
 328   {
 329     return par[0]*(exp(-0.5*t*t));
 330   }
 331
 332   if (t < -absAlpha) //left tail
 333   {
 334     Double_t a =  TMath::Power(par[4]/absAlpha,par[4])*exp(-0.5*absAlpha*absAlpha);
 335     Double_t b = par[4]/absAlpha - absAlpha;
 336     return par[0]*(a/TMath::Power(b - t, par[4]));
 337   }
 338
 339   if (t >= absAlpha2) //right tail
 340   {
 341
 342     Double_t c =  TMath::Power(par[6]/absAlpha2,par[6])*exp(-0.5*absAlpha2*absAlpha2);
 343     Double_t d = par[6]/absAlpha2 - absAlpha2;
 344     return par[0]*(c/TMath::Power(d + t, par[6]));
 345   }
 346
 347   return 0. ;
 348 }
 349
 350 //------------------------------------------------------------------------------
 351 Double_t Gaus(Double_t *x, Double_t *par)
 352 {
 353   // gaussian
 354   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]));
 355
 356 }
 357
 358 //------------------------------------------------------------------------------
 359 Double_t Exp(Double_t *x, Double_t *par)
 360 {
 361   // exponential
 362   return par[0]*TMath::Exp(par[1]*x[0]);
 363
 364 }
 365
 366 //------------------------------------------------------------------------------
 367 Double_t Pow(Double_t *x, Double_t *par)
 368 {
 369   // power law
 370   return par[0]*TMath::Power(x[0],par[1]);
 371
 372 }
 373
 374 //------------------------------------------------------------------------------
 375 Double_t fitFunctionVWG(Double_t *x, Double_t *par)
 376 {
 377   if (x[0] > 2.9 && x[0] < 3.3) TF1::RejectPoint();
 378   return BackgroundVWG(x, par);
 379 }
 380
 381 //------------------------------------------------------------------------------
 382 Double_t fitFunctionCB2VWG(Double_t *x, Double_t *par)
 383 {
 384   return BackgroundVWG(x, par) + CrystalBallExtended(x, &par[4]);
 385 }
 386
 387 //------------------------------------------------------------------------------
 388 Double_t func2CB2VWG(Double_t *x, Double_t *par)
 389 {
 390   /// 2 extended crystal balls + variable width gaussian
 391   /// width of the second CB related to the first (free) one.
 392
 393   Double_t par2[7] = {
 394     par[11],
 395     3.68609+(par[5]-3.096916)/3.096916*3.68609,
 396     par[6]/3.096916*3.68609,
 397     par[7],
 398     par[8],
 399     par[9],
 400     par[10]
 401   };
 402   return BackgroundVWG(x, par) + CrystalBallExtended(x, &par[4]) + CrystalBallExtended(x, par2);
 403 }
 404
 405 //_____________________________________________________________________________
 406 //_____________________________________________________________________________
 407 //_____________________________________________________________________________
 408 //_____________________________________________________________________________
 409 //_____________________________________________________________________________
 410
 411 //_____________________________________________________________________________
 412 AliAnalysisMuMuResult::AliAnalysisMuMuResult(TRootIOCtor* /*io*/) :
 413 TNamed("",""),
 414 fNofRuns(),
 415 fNofTriggers(-1),
 416 fMinv(0x0),
 417 fBin(),
 418 fSubResults(0x0),
 419 fMap(0x0),
 420 fMother(0x0),
 421 fKeys(0x0),
 422 fWeight(0.0),
 423 fRebin(0),
 424 fTriggerClass(),
 425 fEventSelection(),
 426 fPairSelection(),
 427 fCentralitySelection(),
 428 fAlias()
 429 {
 430 }
 431
 432 //_____________________________________________________________________________
 433 AliAnalysisMuMuResult::AliAnalysisMuMuResult(const TH1& hminv)
 434 :
 435   TNamed("",""),
 436   fNofRuns(1),
 437   fNofTriggers(-1),
 438   fMinv(0x0),
 439   fBin(),
 440   fSubResults(0x0),
 441   fMap(0x0),
 442   fMother(0x0),
 443 fKeys(0x0),
 444 fWeight(0.0),
 445 fRebin(0),
 446 fTriggerClass(),
 447 fEventSelection(),
 448 fPairSelection(),
 449 fCentralitySelection(),
 450 fAlias()
 451 {
 452   SetMinv(hminv);
 453 }
 454
 455 //_____________________________________________________________________________
 456 AliAnalysisMuMuResult::AliAnalysisMuMuResult(const TH1& hminv,
 457                                              const char* fitType,
 458                                              Int_t nrebin)
 459 :
 460 TNamed(Form("%s:%d",fitType,nrebin),""),
 461 fNofRuns(1),
 462 fNofTriggers(-1),
 463 fMinv(0x0),
 464 fBin(),
 465 fSubResults(0x0),
 466 fMap(0x0),
 467 fMother(0x0),
 468 fKeys(0x0),
 469 fWeight(0.0),
 470 fRebin(nrebin),
 471 fTriggerClass(),
 472 fEventSelection(),
 473 fPairSelection(),
 474 fCentralitySelection(),
 475 fAlias()
 476 {
 477   SetMinv(hminv);
 478 }
 479
 480 //_____________________________________________________________________________
 481 AliAnalysisMuMuResult::AliAnalysisMuMuResult(const TH1& hminv,
 482                                              const char* triggerName,
 483                                              const char* eventSelection,
 484                                              const char* pairSelection,
 485                                              const char* centSelection,
 486                                              const AliAnalysisMuMuBinning::Range& bin)
 487 :
 488 TNamed(Form("%s-%s-%s-%s",triggerName,eventSelection,pairSelection,centSelection),""),
 489 fNofRuns(1),
 490 fNofTriggers(-1),
 491 fMinv(0x0),
 492 fBin(bin),
 493 fSubResults(0x0),
 494 fMap(0x0),
 495 fMother(0x0),
 496 fKeys(0x0),
 497 fWeight(0.0),
 498 fRebin(1),
 499 fTriggerClass(triggerName),
 500 fEventSelection(eventSelection),
 501 fPairSelection(pairSelection),
 502 fCentralitySelection(centSelection),
 503 fAlias()
 504 {
 505   SetMinv(hminv);
 506 }
 507
 508 //_____________________________________________________________________________
 509 AliAnalysisMuMuResult::AliAnalysisMuMuResult(const AliAnalysisMuMuResult& rhs)
 510 :
 511 TNamed(rhs),
 512 fNofRuns(rhs.NofRuns()),
 513 fNofTriggers(rhs.NofTriggers()),
 514 fMinv(0x0),
 515 fBin(rhs.Bin()),
 516 fSubResults(0x0),
 517 fMap(0x0),
 518 fMother(0x0),
 519 fKeys(0x0),
 520 fWeight(rhs.fWeight),
 521 fRebin(rhs.fRebin),
 522 fAlias()
 523 {
 524   /// copy ctor
 525   /// Note that the mother is lost
 526   /// fKeys remains 0x0 so it will be recomputed if need be
 527
 528   if ( rhs.fMinv )
 529   {
 530     fMinv = static_cast<TH1*>(rhs.fMinv->Clone());
 531   }
 532
 533   if (rhs.fSubResults)
 534   {
 535     fSubResults = static_cast<TObjArray*>(rhs.fSubResults->Clone());
 536   }
 537
 538   if ( rhs.fMap )
 539   {
 540     fMap = static_cast<TMap*>(rhs.fMap->Clone());
 541   }
 542
 543   if ( rhs.fAlias.Length() > 0 )
 544   {
 545     fAlias = rhs.fAlias;
 546   }
 547 }
 548
 549 //_____________________________________________________________________________
 550 AliAnalysisMuMuResult& AliAnalysisMuMuResult::operator=(const AliAnalysisMuMuResult& rhs)
 551 {
 552   /// Assignment operator
 553
 554   if (this!=&rhs)
 555   {
 556     delete fMinv;
 557     delete fMap;
 558     delete fSubResults;
 559
 560     fMinv = 0x0;
 561     fMap = 0x0;
 562     fSubResults = 0x0;
 563     fKeys = 0x0;
 564
 565     if ( rhs.fMinv )
 566     {
 567       fMinv = static_cast<TH1*>(rhs.fMinv->Clone());
 568     }
 569
 570     if (rhs.fSubResults)
 571     {
 572       fSubResults = static_cast<TObjArray*>(rhs.fSubResults->Clone());
 573     }
 574
 575     if ( rhs.fMap )
 576     {
 577       fMap = static_cast<TMap*>(rhs.fMap->Clone());
 578     }
 579
 580     static_cast<TNamed&>(*this)=rhs;
 581
 582     fNofRuns = rhs.NofRuns();
 583     fNofTriggers = rhs.NofTriggers();
 584     fBin = rhs.Bin();
 585     fWeight = rhs.fWeight;
 586     fRebin = rhs.fRebin;
 587
 588     fAlias="";
 589
 590     if ( rhs.fAlias.Length() > 0 )
 591     {
 592       fAlias = rhs.fAlias;
 593     }
 594   }
 595
 596   return *this;
 597 }
 598
 599 //_____________________________________________________________________________
 600 AliAnalysisMuMuResult::~AliAnalysisMuMuResult()
 601 {
 602   // dtor
 603   delete fMap;
 604   delete fMinv;
 605   delete fSubResults;
 606   delete fKeys;
 607 }
 608
 609 //_____________________________________________________________________________
 610 const AliAnalysisMuMuBinning::Range& AliAnalysisMuMuResult::Bin() const
 611 {
 612   /// Get the bin of this result
 613   if ( !Mother() ) return fBin;
 614   else return Mother()->Bin();
 615 }
 616
 617 //_____________________________________________________________________________
 618 TObject* AliAnalysisMuMuResult::Clone(const char* /*newname*/) const
 619 {
 620   /// Clone this result
 621   return new AliAnalysisMuMuResult(*this);
 622 }
 623
 624 //_____________________________________________________________________________
 625 Bool_t AliAnalysisMuMuResult::Correct(const AliAnalysisMuMuResult& other, const char* particle, const char* subResultName)
 626 {
 627   /// Assuming other has an AccxEff entry, correct this value by AccxEff of other
 628
 629   if ( HasValue(Form("Nof%s",particle)) )
 630   {
 631     if (!other.HasValue(Form("AccEff%s",particle),subResultName))
 632     {
 633       AliError(Form("Cannot correct as I do not find the AccEff%s value (subResultName=%s)!",particle,subResultName));
 634       return kFALSE;
 635     }
 636
 637     Double_t acc = other.GetValue(Form("AccEff%s",particle),subResultName);
 638     Double_t value = 0.0;
 639
 640     if ( acc > 0 ) value = GetValue(Form("Nof%s",particle)) / acc;
 641
 642     Double_t error = ErrorAB( GetValue(Form("Nof%s",particle)),
 643                               GetErrorStat(Form("Nof%s",particle)),
 644                               other.GetValue(Form("AccEff%s",particle),subResultName),
 645                               other.GetErrorStat(Form("AccEff%s",particle),subResultName) );
 646
 647     Set(Form("CorrNof%s",particle),value,error*value);
 648
 649     return kTRUE;
 650   }
 651   else
 652   {
 653     AliError(Form("Result does not have Nof%s : cannot correct it !",particle));
 654   }
 655   return kFALSE;
 656 }
 657
 658 //_____________________________________________________________________________
 659 Double_t AliAnalysisMuMuResult::CountParticle(const TH1& hminv, const char* particle, Double_t sigma)
 660 {
 661   /// Count the number of entries in an invariant mass range
 662
 663   const std::map<std::string,Double_t>& m = MassMap();
 664
 665   std::map<std::string,Double_t>::const_iterator it = m.find(particle);
 666
 667   if ( it == m.end() )
 668   {
 669     AliErrorClass(Form("Don't know the mass of particle %s",particle));
 670     return 0.0;
 671   }
 672
 673   Double_t mass = it->second;
 674
 675   if ( sigma < 0 )
 676   {
 677     AliDebugClass(1,Form("Oups. Got a sigma of %e for particle %s !",sigma,particle));
 678     return hminv.Integral();
 679   }
 680
 681   TAxis* x = hminv.GetXaxis();
 682
 683   Int_t b1 = x->FindBin(mass-sigma);
 684   Int_t b2 = x->FindBin(mass+sigma);
 685
 686   AliDebugClass(1,Form("hminv getentries %e integral %e",hminv.GetEntries(),hminv.Integral(b1,b2)));
 687
 688   return hminv.Integral(b1,b2);
 689 }
 690
 691 /*
 692 //_____________________________________________________________________________
 693 void AliAnalysisMuMuResult::FitJpsiPsiPrimeCustom(TH1& h)
 694 {
 695   std::cout << "Fit with jpsi + psiprime (custom)" << std::endl;
 696
 697   const Double_t xmin(1.5);
 698   const Double_t xmax(8.0);
 699
 700   fitTotal = new TF1("fitTotal",funcJpsiPsiPrimeCustom,xmin,xmax,14);
 701   fitTotal->SetLineColor(4);
 702
 703   fitTotal->SetParName(0,"cstecb");
 704   fitTotal->SetParName(1,"alpha");
 705   fitTotal->SetParName(2,"n");
 706   fitTotal->SetParName(3,"meanjpsi");
 707   fitTotal->SetParName(4,"sigmajpsi");
 708   fitTotal->SetParName(5,"alphaprime");
 709   fitTotal->SetParName(6,"nprime");
 710   fitTotal->SetParName(7,"cstepol1");
 711   fitTotal->SetParName(8,"slopepol1");
 712   fitTotal->SetParName(9,"cstegaus");
 713   fitTotal->SetParName(10,"meanpsiprime");
 714   fitTotal->SetParName(11,"sigmapsiprime");
 715   fitTotal->SetParName(12,"csteexpo");
 716   fitTotal->SetParName(13,"slopeexpo");
 717
 718   fitTotal->SetParameter( 0,1);
 719
 720   const char* fitOption = "SQBR+";
 721   const Double_t alphaMC = 0.936;
 722   const Double_t nMC = 4.44;
 723   const Double_t alphaprimeMC = 1.60;
 724   const Double_t nprimeMC = 3.23;
 725
 726   TF1* fcb = new TF1("cb",funcCB2,2.9,3.3,7);
 727   fcb->SetParameters(1,1.0,4.0,3.1,0.1,1.5,3);
 728
 729   fcb->SetParLimits(3,3,4);
 730   fcb->SetParLimits(4,0,1);
 731
 732   fcb->FixParameter(1,alphaMC);
 733   fcb->FixParameter(2,nMC);
 734   fcb->FixParameter(5,alphaprimeMC);
 735   fcb->FixParameter(6,nprimeMC);
 736
 737   TFitResultPtr rcb = h.Fit(fcb,fitOption,"",2.9,3.3);
 738
 739   if (!rcb.Get())
 740   {
 741     return;
 742   }
 743
 744   fitTotal->SetParameter(0,rcb->Parameter(0));
 745   fitTotal->SetParameter(1,rcb->Parameter(1)); fitTotal->SetParLimits(1,0,10); // alpha
 746   fitTotal->SetParameter(2,rcb->Parameter(2)); fitTotal->SetParLimits(2,1,10); // n
 747   fitTotal->SetParameter(3,rcb->Parameter(3)); fitTotal->SetParLimits(3,3.0,3.5); // mean
 748   fitTotal->SetParameter(4,rcb->Parameter(4)); fitTotal->SetParLimits(4,0,1); // sigma
 749   fitTotal->SetParameter(5,rcb->Parameter(5)); fitTotal->SetParLimits(1,0,10); // alphaprime
 750   fitTotal->SetParameter(6,rcb->Parameter(6)); fitTotal->SetParLimits(2,1,10); // nprime
 751
 752   fitTotal->FixParameter(1,alphaMC);
 753   fitTotal->FixParameter(2,nMC);
 754   fitTotal->FixParameter(5,alphaprimeMC);
 755   fitTotal->FixParameter(6,nprimeMC);
 756
 757   TF1* fge = new TF1("fge","gaus(0)+expo(3)",3.5,4.4);
 758   fge->SetParameters(1,3.6,0.25,1,1);
 759   TFitResultPtr rpsiprime = h.Fit(fge,fitOption,"",3.5,4.4);
 760
 761   if (static_cast<int>(rpsiprime))
 762   {
 763     AliInfo("Will fix psiprime parameters");
 764     fitTotal->FixParameter(9,0);
 765     fitTotal->FixParameter(10,3.7);
 766     fitTotal->FixParameter(11,0.1);
 767   }
 768   else
 769   {
 770     fitTotal->SetParameter(10,rpsiprime->Parameter(1)); fitTotal->SetParLimits(10,3.5,3.8); // mean'
 771     fitTotal->SetParameter(11,rpsiprime->Parameter(2)); fitTotal->SetParLimits(11,0.05,0.7); // sigma'
 772   }
 773
 774   TFitResultPtr rpol1 = h.Fit("pol1",fitOption,"",1.5,2.5);
 775   fitTotal->SetParameter( 7,rpol1->Parameter(0));
 776   fitTotal->SetParameter( 8,rpol1->Parameter(1));
 777
 778   TFitResultPtr rexpo = h.Fit("expo",fitOption,"",4.5,7.0);
 779   fitTotal->SetParameter(12,rexpo->Parameter(0));
 780   fitTotal->SetParameter(13,rexpo->Parameter(1));
 781
 782
 783   TFitResultPtr r = h.Fit(fitTotal,fitOption,"",1.5,7);
 784
 785   TF1* signal = new TF1("signal","gaus",2,6);
 786   signal->SetParameters(fitTotal->GetParameter(0),
 787                         fitTotal->GetParameter(3),
 788                         fitTotal->GetParameter(4));
 789
 790   TF1* signalPrime = new TF1("signalPrime","gaus",2,6);
 791   signalPrime->SetParameters(fitTotal->GetParameter(9),
 792                              fitTotal->GetParameter(10),
 793                              fitTotal->GetParameter(11));
 794
 795   Double_t gausParameters[3];
 796   Double_t covarianceMatrix[3][3];
 797   Double_t gausParametersPrime[3];
 798   Double_t covarianceMatrixPrime[3][3];
 799
 800   covarianceMatrix[0][0] = (r->GetCovarianceMatrix())(0,0);
 801   covarianceMatrix[1][0] = (r->GetCovarianceMatrix())(3,0);
 802   covarianceMatrix[2][0] = (r->GetCovarianceMatrix())(4,0);
 803   covarianceMatrix[0][1] = (r->GetCovarianceMatrix())(0,3);
 804   covarianceMatrix[0][2] = (r->GetCovarianceMatrix())(0,4);
 805
 806   for ( int iy = 1; iy < 3; ++iy )
 807   {
 808     for ( int ix = 1; ix < 3; ++ix )
 809     {
 810       covarianceMatrix[ix][iy] = (r->GetCovarianceMatrix())(ix+2,iy+2);
 811     }
 812   }
 813
 814   gausParameters[0] = fitTotal->GetParameter(0);
 815   gausParameters[1] = fitTotal->GetParameter(3);
 816   gausParameters[2] = fitTotal->GetParameter(4);
 817
 818   gausParametersPrime[0] = fitTotal->GetParameter(9);
 819   gausParametersPrime[1] = fitTotal->GetParameter(10);
 820   gausParametersPrime[2] = fitTotal->GetParameter(11);
 821
 822   covarianceMatrixPrime[0][0] = (r->GetCovarianceMatrix())(9,9);
 823   covarianceMatrixPrime[1][0] = (r->GetCovarianceMatrix())(10,9);
 824   covarianceMatrixPrime[2][0] = (r->GetCovarianceMatrix())(11,9);
 825   covarianceMatrixPrime[0][1] = (r->GetCovarianceMatrix())(9,10);
 826   covarianceMatrixPrime[0][2] = (r->GetCovarianceMatrix())(9,11);
 827
 828   for ( int iy = 1; iy < 3; ++iy )
 829   {
 830     for ( int ix = 1; ix < 3; ++ix )
 831     {
 832       covarianceMatrixPrime[ix][iy] = (r->GetCovarianceMatrix())(ix+2,iy+2);
 833     }
 834   }
 835
 836   double n = signal->Integral(2,6)/h.GetBinWidth(10);
 837   double nerr = signal->IntegralError(2,6,&gausParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(10);
 838   Set("NofJpsi",n,nerr);
 839   Set("MeanJpsi",fitTotal->GetParameter(3),fitTotal->GetParError(3));
 840   Set("SigmaJpsi",fitTotal->GetParameter(4),fitTotal->GetParError(4));
 841
 842   double nprime = signalPrime->Integral(2,6)/h.GetBinWidth(10);
 843   double nerrprime = signalPrime->IntegralError(2,6,&gausParametersPrime[0],&covarianceMatrixPrime[0][0])/h.GetBinWidth(10);
 844   Set("NofPsiPrime",nprime,nerrprime);
 845   Set("MeanPsiPrime",fitTotal->GetParameter(10),fitTotal->GetParError(10));
 846   Set("SigmaPsiPrime",fitTotal->GetParameter(11),fitTotal->GetParError(11));
 847 }
 848 */
 849
 850 /*
 851 //_____________________________________________________________________________
 852 AliAnalysisMuMuResult::SubResult AliAnalysisMuMuResult::FitJpsiPsiPrimeCB(TH1& h)
 853 {
 854   std::cout << "Fit with jpsi + psiprime (CB) " << std::endl;
 855
 856   const Double_t xmin(1.5);
 857   const Double_t xmax(8.0);
 858
 859   fitTotal = new TF1("fitTotal",funcJpsiPsiPrime,xmin,xmax,14);
 860
 861 //  Double_t N = par[0];
 862 //  Double_t alpha = par[1];
 863 //  Double_t n = par[2];
 864 //  Double_t mean = par[3];
 865 //  Double_t sigma = par[4];
 866
 867   fitTotal->SetParameter( 0,1); // N
 868   fitTotal->FixParameter( 1,0.936); // alpha
 869   fitTotal->FixParameter( 2,4.44); // n
 870   fitTotal->SetParameter( 3,3.1); fitTotal->SetParLimits(3,3.0,3.2); // mean
 871   fitTotal->SetParameter( 4,0.07); fitTotal->SetParLimits(4,0.02,1); // sigma
 872
 873   fitTotal->SetParameter( 5,0.01); // N'
 874   fitTotal->FixParameter( 6,0.936); // alpha'
 875   fitTotal->FixParameter( 7,4.44); // n'
 876   fitTotal->SetParameter( 8,3.7); fitTotal->SetParLimits(8,3.5,3.8); // mean'
 877   fitTotal->SetParameter( 9,0.1); fitTotal->SetParLimits(9,0.02,1.0); // sigma'
 878
 879   fitTotal->SetParameter(10,h.GetMaximum());
 880   fitTotal->SetParameter(11,-1);
 881
 882   fitTotal->SetParameter(12,h.GetMaximum()/100);
 883   fitTotal->SetParameter(13,-1);
 884
 885   TFitResultPtr r = h.Fit(fitTotal,"SQBI","",1.5,6);
 886
 887 //  for ( int ix = 0; ix < fitTotal->GetNpar(); ++ix )
 888 //  {
 889 //    for ( int iy = 0; iy < fitTotal->GetNpar(); ++iy )
 890 //    {
 891 //      std::cout << Form("COV(%d,%d)=%e ",ix,iy,r->GetCovarianceMatrix()(ix,iy));
 892 //    }
 893 //    std::cout << std::endl;
 894 //  }
 895
 896
 897   TF1* signal = new TF1("signal","gaus",2,8);
 898
 899   signal->SetParameters(fitTotal->GetParameter(0),
 900                         fitTotal->GetParameter(3),
 901                         fitTotal->GetParameter(4));
 902
 903   TF1* signalPrime = new TF1("signalPrime","gaus",2,8);
 904
 905   signalPrime->SetParameters(fitTotal->GetParameter(0),
 906                              fitTotal->GetParameter(8),
 907                              fitTotal->GetParameter(9));
 908
 909   Double_t gausParameters[3];
 910   Double_t gausParametersPrime[3];
 911   Double_t covarianceMatrix[3][3];
 912   Double_t covarianceMatrixPrime[3][3];
 913
 914   gausParameters[0] = fitTotal->GetParameter(0);
 915   gausParameters[1] = fitTotal->GetParameter(3);
 916   gausParameters[2] = fitTotal->GetParameter(4);
 917
 918   covarianceMatrix[0][0] = (r->GetCovarianceMatrix())(0,0);
 919   covarianceMatrix[1][0] = (r->GetCovarianceMatrix())(3,0);
 920   covarianceMatrix[2][0] = (r->GetCovarianceMatrix())(4,0);
 921   covarianceMatrix[0][1] = (r->GetCovarianceMatrix())(0,3);
 922   covarianceMatrix[0][2] = (r->GetCovarianceMatrix())(0,4);
 923
 924   for ( int iy = 1; iy < 3; ++iy )
 925   {
 926     for ( int ix = 1; ix < 3; ++ix )
 927     {
 928       covarianceMatrix[ix][iy] = (r->GetCovarianceMatrix())(ix+2,iy+2);
 929     }
 930   }
 931
 932   gausParametersPrime[0] = fitTotal->GetParameter(0);
 933   gausParametersPrime[1] = fitTotal->GetParameter(8);
 934   gausParametersPrime[2] = fitTotal->GetParameter(9);
 935
 936   covarianceMatrixPrime[0][0] = (r->GetCovarianceMatrix())(0,0);
 937   covarianceMatrixPrime[1][0] = (r->GetCovarianceMatrix())(8,0);
 938   covarianceMatrixPrime[2][0] = (r->GetCovarianceMatrix())(9,0);
 939   covarianceMatrixPrime[0][1] = (r->GetCovarianceMatrix())(0,8);
 940   covarianceMatrixPrime[0][2] = (r->GetCovarianceMatrix())(0,9);
 941
 942   for ( int iy = 1; iy < 3; ++iy )
 943   {
 944     for ( int ix = 1; ix < 3; ++ix )
 945     {
 946       covarianceMatrixPrime[ix][iy] = (r->GetCovarianceMatrix())(ix+2,iy+2);
 947     }
 948   }
 949
 950   double n = signal->Integral(2,6)/h.GetBinWidth(10);
 951   double nerr = signal->IntegralError(2,6,&gausParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(10);
 952
 953   Set("NofJpsi",n,nerr);
 954   Set("MeanJpsi",fitTotal->GetParameter(3),fitTotal->GetParError(3));
 955   Set("SigmaJpsi",fitTotal->GetParameter(4),fitTotal->GetParError(4));
 956
 957   double nprime = signalPrime->Integral(2,6)/h.GetBinWidth(10);
 958   double nerrprime = signalPrime->IntegralError(2,6,&gausParametersPrime[0],&covarianceMatrixPrime[0][0])/h.GetBinWidth(10);
 959
 960   Set("NofPsiPrime",nprime,nerrprime);
 961   Set("MeanPsiPrime",fitTotal->GetParameter(8),fitTotal->GetParError(8));
 962   Set("SigmaPsiPrime",fitTotal->GetParameter(9),fitTotal->GetParError(9));
 963
 964 }
 965   */
 966
 967 //_____________________________________________________________________________
 968 AliAnalysisMuMuResult* AliAnalysisMuMuResult::FitJpsiGCBE(TH1& h)
 969 {
 970   /// Fit Jpsi spectra with crystal ball + gaussian + exponential
 971
 972   std::cout << "Fit with jpsi alone (gaus + CB + expo)" << std::endl;
 973
 974   Int_t nrebin = fMinv->GetXaxis()->GetNbins() / h.GetXaxis()->GetNbins();
 975
 976   AliAnalysisMuMuResult* r = new AliAnalysisMuMuResult(h,"JPSIGCBE",nrebin);
 977
 978   TH1* hfit = r->Minv();
 979
 980   const Double_t xmin(1.0);
 981   const Double_t xmax(8.0);
 982
 983   TF1* fitTotal = new TF1("fitTotal",funcJpsiGCBE,xmin,xmax,10);
 984   fitTotal->SetParNames("cste","x0","sigma0","N","alpha","n","mean","sigma","expocste","exposlope");
 985
 986   fitTotal->SetParLimits(3,0,h.GetMaximum()*2); // N
 987
 988   const Double_t cbalpha(0.98);
 989   const Double_t cbn(5.2);
 990
 991   fitTotal->FixParameter(4,cbalpha);
 992   fitTotal->FixParameter(5,cbn);
 993
 994   fitTotal->SetParLimits(6,2.8,3.2); // mean
 995   fitTotal->SetParLimits(7,0.02,0.3); // sigma
 996
 997   TF1* fg = new TF1("fg","gaus",xmin,xmax);
 998
 999   hfit->Fit(fg,"","",0.75,3.0);
1000
1001   fitTotal->SetParameter(0,fg->GetParameter(0));
1002   fitTotal->SetParameter(1,fg->GetParameter(1));
1003   fitTotal->SetParameter(2,fg->GetParameter(2));
1004
1005   TF1* fexpo = new TF1("expo","expo",xmin,xmax);
1006
1007   hfit->Fit(fexpo,"","",3.5,5);
1008
1009   fitTotal->SetParameter(8,fexpo->GetParameter(0));
1010   fitTotal->SetParameter(9,fexpo->GetParameter(1));
1011
1012   fitTotal->SetParameter(3,h.GetMaximum()),
1013   fitTotal->SetParameter(4,cbalpha);
1014   fitTotal->SetParameter(5,cbn);
1015   fitTotal->SetParameter(6,3.15);
1016   fitTotal->SetParameter(7,0.1);
1017
1018   const char* fitOption = "SI+";
1019
1020   TFitResultPtr fitResult = hfit->Fit(fitTotal,fitOption,"",2,5);
1021
1022   r->Set("MeanJpsi",fitTotal->GetParameter(6),fitTotal->GetParError(6));
1023   r->Set("SigmaJpsi",fitTotal->GetParameter(7),fitTotal->GetParError(7));
1024
1025   double m = r->GetValue("MeanJpsi");
1026   double s = r->GetValue("SigmaJpsi");
1027   double n = 3.0;
1028
1029   TF1* fcb = new TF1("fcb",funcCB,xmin,xmax,5);
1030   fcb->SetParameters(fitTotal->GetParameter(3),
1031                      fitTotal->GetParameter(4),
1032                      fitTotal->GetParameter(5),
1033                      fitTotal->GetParameter(6),
1034                      fitTotal->GetParameter(7));
1035
1036   fcb->SetLineColor(6);
1037   fcb->SetNpx(100);
1038   TLine* l1 = new TLine(m-n*s,0,m-n*s,fitTotal->GetParameter(3));
1039   TLine* l2 = new TLine(m+n*s,0,m+n*s,fitTotal->GetParameter(3));
1040   l1->SetLineColor(6);
1041   l2->SetLineColor(6);
1042   h.GetListOfFunctions()->Add(fcb);
1043   h.GetListOfFunctions()->Add(l1);
1044   h.GetListOfFunctions()->Add(l2);
1045
1046
1047   Double_t cbParameters[5];
1048   Double_t covarianceMatrix[5][5];
1049
1050   cbParameters[0] = fitTotal->GetParameter(3);
1051   cbParameters[1] = fitTotal->GetParameter(4);
1052   cbParameters[2] = fitTotal->GetParameter(5);
1053   cbParameters[3] = fitTotal->GetParameter(6);
1054   cbParameters[4] = fitTotal->GetParameter(7);
1055
1056   for ( int iy = 0; iy < 5; ++iy )
1057   {
1058     for ( int ix = 0; ix < 5; ++ix )
1059     {
1060       covarianceMatrix[ix][iy] = (fitResult->GetCovarianceMatrix())(ix+3,iy+3);
1061     }
1062   }
1063
1064   double njpsi = fcb->Integral(m-n*s,m+n*s)/h.GetBinWidth(1);
1065
1066   double nerr = fcb->IntegralError(m-n*s,m+n*s,&cbParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(1);
1067
1068   r->Set("NofJpsi",njpsi,nerr);
1069
1070   return r;
1071 }
1072
1073 /*
1074 //_____________________________________________________________________________
1075 void AliAnalysisMuMuResult::FitJpsiPCBE(TH1& h)
1076 {
1077   std::cout << "Fit with jpsi alone (pol2 + CB + expo)" << std::endl;
1078
1079   const Double_t xmin(2.0);
1080   const Double_t xmax(5.0);
1081
1082   fitTotal = new TF1("fitTotal",funcJpsiPCBE,xmin,xmax,10);
1083   fitTotal->SetParNames("p0","p1","p2","N","alpha","n","mean","sigma","expocste","exposlope");
1084
1085   fitTotal->SetParLimits(3,0,h.GetMaximum()*2); // N
1086
1087   const Double_t cbalpha(0.98);
1088   const Double_t cbn(5.2);
1089
1090   fitTotal->FixParameter(4,cbalpha);
1091   fitTotal->FixParameter(5,cbn);
1092
1093   fitTotal->SetParLimits(6,2,4); // mean
1094   fitTotal->SetParLimits(7,0.05,0.2); // sigma
1095
1096   TF1* fpol2 = new TF1("pol2","pol2",xmin,xmax);
1097
1098   h.Fit(fpol2,"+","",2,2.8);
1099
1100   fitTotal->SetParameter(0,fpol2->GetParameter(0));
1101   fitTotal->SetParameter(1,fpol2->GetParameter(1));
1102   fitTotal->SetParameter(2,fpol2->GetParameter(2));
1103
1104   TF1* fexpo = new TF1("expo","expo",xmin,xmax);
1105
1106   h.Fit(fexpo,"+","",3.5,4.5);
1107
1108   fitTotal->SetParameter(8,fexpo->GetParameter(0));
1109   fitTotal->SetParameter(9,fexpo->GetParameter(1));
1110
1111   fitTotal->SetParameter(3,h.GetMaximum()),
1112   fitTotal->SetParameter(4,cbalpha);
1113   fitTotal->SetParameter(5,cbn);
1114   fitTotal->SetParameter(6,3.15);
1115   fitTotal->SetParameter(7,0.1);
1116
1117   h.Fit(fitTotal,"+","",2.5,5);
1118
1119   Set("MeanJpsi",fitTotal->GetParameter(6),fitTotal->GetParError(6));
1120   Set("SigmaJpsi",fitTotal->GetParameter(7),fitTotal->GetParError(7));
1121
1122   double m = GetValue("MeanJpsi");
1123   double s = GetValue("SigmaJpsi");
1124   double n = 2.0;
1125
1126   TF1* fcb = new TF1("fcb",funcCB,xmin,xmax,5);
1127   fcb->SetParameters(fitTotal->GetParameter(3),
1128                      fitTotal->GetParameter(4),
1129                      fitTotal->GetParameter(5),
1130                      fitTotal->GetParameter(6),
1131                      fitTotal->GetParameter(7));
1132
1133   fcb->SetLineColor(6);
1134   fcb->SetNpx(100);
1135   TLine* l1 = new TLine(m-n*s,0,m-n*s,fitTotal->GetParameter(3));
1136   TLine* l2 = new TLine(m+n*s,0,m+n*s,fitTotal->GetParameter(3));
1137   l1->SetLineColor(6);
1138   l2->SetLineColor(6);
1139   h.GetListOfFunctions()->Add(fcb);
1140   h.GetListOfFunctions()->Add(l1);
1141   h.GetListOfFunctions()->Add(l2);
1142
1143
1144   Set("NofJpsi",fitTotal->Integral(m-n*s,m+n*s)/h.GetBinWidth(1),fitTotal->IntegralError(m-n*s,m+n*s)/h.GetBinWidth(1));
1145
1146   //  Set("NofJpsi",fitTotal->Integral(0,10)/h.GetBinWidth(1),fitTotal->IntegralError(0,10)/h.GetBinWidth(1));
1147
1148 }
1149
1150 //_____________________________________________________________________________
1151 void AliAnalysisMuMuResult::FitJpsiCBE(TH1& h)
1152 {
1153   std::cout << "Fit with jpsi alone" << std::endl;
1154
1155   const Double_t xmin(1.5);
1156   const Double_t xmax(8.0);
1157
1158   fitTotal = new TF1("fitTotal",funcJpsiCBE,xmin,xmax,7);
1159   fitTotal->SetParNames("N","alpha","n","mean","sigma","expocste","exposlope");
1160
1161 //  fitTotal->SetParameters(h.GetMaximum(),1,5,3.0,0.07,1.5,3,1,0);
1162
1163   fitTotal->SetParameters(1,1.15,3.6,3.0,0.07,1,-1);
1164
1165   fitTotal->SetParLimits(0,0,h.GetMaximum()); // N
1166 //  fitTotal->SetParLimits(1,0.1,2); // alpha
1167   fitTotal->FixParameter(1,0.98);
1168 //  fitTotal->SetParLimits(2,0.01,5); // n
1169   fitTotal->FixParameter(2,5.2);
1170   fitTotal->SetParLimits(3,2.8,3.5); // mean
1171   fitTotal->SetParLimits(4,0.05,0.2); // sigma
1172
1173   TF1* fexpo = new TF1("expo","expo",xmin,xmax);
1174
1175   h.Fit(fexpo,"+","",2,3);
1176
1177   fitTotal->SetParameter(5,fexpo->GetParameter(0));
1178   fitTotal->SetParameter(6,fexpo->GetParameter(1));
1179
1180   h.Fit(fitTotal,"+","",2,5);
1181
1182
1183   Set("MeanJpsi",fitTotal->GetParameter(3),fitTotal->GetParError(3));
1184   Set("SigmaJpsi",fitTotal->GetParameter(4),fitTotal->GetParError(4));
1185
1186   double m = GetValue("MeanJpsi");
1187   double s = GetValue("SigmaJpsi");
1188   double n = 3.0;
1189
1190   TF1* fcb = new TF1("fcb",funcCB,xmin,xmax,5);
1191   fcb->SetParameters(fitTotal->GetParameter(0),
1192                      fitTotal->GetParameter(1),
1193                      fitTotal->GetParameter(2),
1194                      fitTotal->GetParameter(3),
1195                      fitTotal->GetParameter(4));
1196
1197   fcb->SetLineColor(6);
1198   fcb->SetNpx(1000);
1199   TLine* l1 = new TLine(m-n*s,0,m-n*s,fitTotal->GetParameter(0));
1200   TLine* l2 = new TLine(m+n*s,0,m+n*s,fitTotal->GetParameter(0));
1201   l1->SetLineColor(6);
1202   l2->SetLineColor(6);
1203   h.GetListOfFunctions()->Add(fcb);
1204   h.GetListOfFunctions()->Add(l1);
1205   h.GetListOfFunctions()->Add(l2);
1206
1207
1208   Set("NofJpsi",fitTotal->Integral(m-n*s,m+n*s)/h.GetBinWidth(1),fitTotal->IntegralError(m-n*s,m+n*s)/h.GetBinWidth(1));
1209
1210   //  Set("NofJpsi",fitTotal->Integral(0,10)/h.GetBinWidth(1),fitTotal->IntegralError(0,10)/h.GetBinWidth(1));
1211
1212 }
1213
1214 //_____________________________________________________________________________
1215 void AliAnalysisMuMuResult::FitJpsiECBE(TH1& h)
1216 {
1217   std::cout << "Fit with jpsi alone (expo + CB + expo)" << std::endl;
1218
1219   const Double_t xmin(1.5);
1220   const Double_t xmax(8.0);
1221
1222   fitTotal = new TF1("fitTotal",funcJpsiECBE,xmin,xmax,9);
1223   fitTotal->SetParNames("e0","s0","N","alpha","n","mean","sigma","e1","s1");
1224
1225   fitTotal->SetParameters(1,-1,1,1.15,3.6,3.2,0.06,-1);
1226
1227   fitTotal->SetParLimits(0,0,h.GetMaximum()*2);
1228
1229   fitTotal->FixParameter(3,0.98); // alpha
1230   fitTotal->FixParameter(4,5.2); // n
1231   fitTotal->SetParLimits(5,2.8,3.5); // mean
1232   fitTotal->SetParLimits(6,0.05,0.2); // sigma
1233
1234   TF1* fexpo1 = new TF1("expo1","expo",xmin,xmax);
1235   TF1* fexpo2 = new TF1("expo2","expo",xmin,xmax);
1236
1237   h.Fit(fexpo1,"","",1.5,3);
1238
1239   fitTotal->SetParameter(0,fexpo1->GetParameter(0));
1240   fitTotal->SetParameter(1,fexpo1->GetParameter(1));
1241
1242   h.Fit(fexpo2,"","",3.5,5.0);
1243
1244   fitTotal->SetParameter(7,fexpo2->GetParameter(0));
1245   fitTotal->SetParameter(8,fexpo2->GetParameter(1));
1246
1247   const char* fitOption = "SI+";
1248
1249   TFitResultPtr r = h.Fit(fitTotal,fitOption,"",2,5);
1250
1251   Set("MeanJpsi",fitTotal->GetParameter(5),fitTotal->GetParError(5));
1252   Set("SigmaJpsi",fitTotal->GetParameter(6),fitTotal->GetParError(6));
1253
1254   double m = GetValue("MeanJpsi");
1255   double s = GetValue("SigmaJpsi");
1256   double n = 3.0;
1257
1258   TF1* fcb = new TF1("fcb",funcCB,xmin,xmax,5);
1259   fcb->SetParameters(fitTotal->GetParameter(2),
1260                      fitTotal->GetParameter(3),
1261                      fitTotal->GetParameter(4),
1262                      fitTotal->GetParameter(5),
1263                      fitTotal->GetParameter(6));
1264
1265   fcb->SetParError(0,fitTotal->GetParError(2));
1266   fcb->SetParError(1,fitTotal->GetParError(3));
1267   fcb->SetParError(2,fitTotal->GetParError(4));
1268   fcb->SetParError(3,fitTotal->GetParError(5));
1269   fcb->SetParError(4,fitTotal->GetParError(6));
1270
1271   fcb->SetLineColor(6);
1272   fcb->SetNpx(1000);
1273   TLine* l1 = new TLine(m-n*s,0,m-n*s,fitTotal->GetParameter(2));
1274   TLine* l2 = new TLine(m+n*s,0,m+n*s,fitTotal->GetParameter(2));
1275   l1->SetLineColor(6);
1276   l2->SetLineColor(6);
1277   h.GetListOfFunctions()->Add(fcb);
1278   h.GetListOfFunctions()->Add(l1);
1279   h.GetListOfFunctions()->Add(l2);
1280
1281
1282   Double_t cbParameters[5];
1283   Double_t covarianceMatrix[5][5];
1284
1285   cbParameters[0] = fitTotal->GetParameter(2);
1286   cbParameters[1] = fitTotal->GetParameter(3);
1287   cbParameters[2] = fitTotal->GetParameter(4);
1288   cbParameters[3] = fitTotal->GetParameter(5);
1289   cbParameters[4] = fitTotal->GetParameter(6);
1290
1291   for ( int iy = 0; iy < 5; ++iy )
1292   {
1293     for ( int ix = 0; ix < 5; ++ix )
1294     {
1295       covarianceMatrix[ix][iy] = (r->GetCovarianceMatrix())(ix+2,iy+2);
1296     }
1297   }
1298
1299
1300   double njpsi = fcb->Integral(m-n*s,m+n*s)/h.GetBinWidth(1);
1301
1302   double nerr = fcb->IntegralError(m-n*s,m+n*s,&cbParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(1);
1303
1304
1305   Set("NofJpsi",njpsi,nerr);
1306 }
1307
1308  */
1309
1310 //_____________________________________________________________________________
1311 AliAnalysisMuMuResult* AliAnalysisMuMuResult::FitJpsi(TH1& h)
1312 {
1313   /// Fit Jpsi spectra using extended crystall ball (CB2) with free tails
1314
1315   std::cout << "Fit with jpsi alone" << std::endl;
1316
1317   Int_t nrebin = fMinv->GetXaxis()->GetNbins() / h.GetXaxis()->GetNbins();
1318
1319   AliAnalysisMuMuResult* r = new AliAnalysisMuMuResult(h,"JPSI",nrebin);
1320
1321   TH1* hfit = r->Minv();
1322
1323   const Double_t xmin(1.5);
1324   const Double_t xmax(8.0);
1325
1326   TF1* fitTotal = new TF1("fitTotal",funcCB2,xmin,xmax,7);
1327   fitTotal->SetParNames("N","alphaLow","nLow","mean","sigma","alphaUp","nUp");
1328   fitTotal->SetParameters(h.GetMaximum(),1,5,3.1,0.07,1.5,3);
1329   fitTotal->SetParLimits(0,0,h.GetMaximum()*2); // N
1330   fitTotal->SetParLimits(1,0,10); // alpha
1331   fitTotal->SetParLimits(2,0.1,10); // n
1332   fitTotal->SetParLimits(3,3,3.1); // mean
1333   fitTotal->SetParLimits(4,0.01,1); // sigma
1334   fitTotal->SetParLimits(5,0,10); // alpha
1335   fitTotal->SetParLimits(6,0.1,10); // n
1336
1337   hfit->Fit(fitTotal,"+","",2,5);
1338
1339
1340   r->Set("MeanJpsi",fitTotal->GetParameter(3),fitTotal->GetParError(3));
1341   r->Set("SigmaJpsi",fitTotal->GetParameter(4),fitTotal->GetParError(4));
1342
1343   double m = r->GetValue("MeanJpsi");
1344   double s = r->GetValue("SigmaJpsi");
1345   double n = 10.0;
1346
1347   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));
1348
1349   return r;
1350 }
1351
1352 //_____________________________________________________________________________
1353 AliAnalysisMuMuResult* AliAnalysisMuMuResult::FitJpsiCB2VWG(const TH1& h)
1354 {
1355   /// Fit Jpsi spectra using extended crystal ball (CB2) + variable width gaussian (VWG)
1356
1357   std::cout << "Fit with jpsi VWG" << std::endl;
1358
1359   Int_t nrebin = fMinv->GetXaxis()->GetNbins() / h.GetXaxis()->GetNbins();
1360
1361   AliAnalysisMuMuResult* r = new AliAnalysisMuMuResult(h,"JPSICB2VWG",nrebin);
1362
1363
1364   TH1* hfit = r->Minv();
1365
1366   const Double_t xmin(2.0);
1367   const Double_t xmax(5.0);
1368
1369 //  // gaussian variable width
1370 //  Double_t sigma = par[2]+par[3]*((x[0]-par[1])/par[1]);
1371 //  return par[0]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/(2.*sigma*sigma));
1372 //  Double_t CrystalBallExtended(Double_t *x,Double_t *par)
1373 //  //par[0] = Normalization
1374 //  //par[1] = mean
1375 //  //par[2] = sigma
1376 //  //par[3] = alpha
1377 //  //par[4] = n
1378 //  //par[5] = alpha'
1379 //  //par[6] = n'
1380
1381   TF1* fitTotal = new TF1("fitTotal",fitFunctionCB2VWG,xmin,xmax,11);
1382   fitTotal->SetParNames("kVWG","mVWG","sVWG1","sVWG2","norm","mean","sigma","alpha","n","alpha'","n'");
1383
1384   fitTotal->SetParameter(0, 10000.); // kVWG
1385   fitTotal->SetParameter(1, 1.9); // mVWG
1386
1387   fitTotal->SetParameter(2, 0.5); // sVWG1
1388   fitTotal->SetParLimits(2, 0., 100.);
1389
1390   fitTotal->SetParameter(3, 0.3); // sVWG2
1391   fitTotal->SetParLimits(3, 0., 100.);
1392
1393   fitTotal->SetParameter(4, h.GetMaximum()); // norm
1394
1395   fitTotal->SetParameter(5, 3.1); // mean
1396   fitTotal->SetParLimits(5, 3.0, 3.2);
1397
1398   fitTotal->SetParameter(6, 0.08); // sigma
1399   fitTotal->SetParLimits(6, 0.04, 0.20);
1400
1401   fitTotal->SetParameter(7,1.0); // alpha
1402   fitTotal->SetParameter(8,5); // n
1403   fitTotal->SetParameter(9,2.0); // alpha'
1404   fitTotal->SetParameter(10,4); // n'
1405
1406 //  fitTotal->FixParameter(7, 0.93);
1407 //  fitTotal->FixParameter(8, 5.59);
1408 //  fitTotal->FixParameter(9, 2.32);
1409 //  fitTotal->FixParameter(10, 3.39);
1410 //  fitTotal->SetParameter(11, 10.);
1411
1412   const char* fitOption = "SIER"; //+";
1413
1414   TFitResultPtr fitResult = hfit->Fit(fitTotal,fitOption,"");
1415
1416   r->Set("MeanJpsi",fitTotal->GetParameter(5),fitTotal->GetParError(5));
1417   r->Set("SigmaJpsi",fitTotal->GetParameter(6),fitTotal->GetParError(6));
1418
1419   double m = r->GetValue("MeanJpsi");
1420   double s = r->GetValue("SigmaJpsi");
1421   double n = 3.0;
1422
1423   TF1* fcb = new TF1("fcb",CrystalBallExtended,xmin,xmax,7);
1424   fcb->SetParameters(fitTotal->GetParameter(4),
1425                      fitTotal->GetParameter(5),
1426                      fitTotal->GetParameter(6),
1427                      fitTotal->GetParameter(7),
1428                      fitTotal->GetParameter(8),
1429                      fitTotal->GetParameter(9),
1430                      fitTotal->GetParameter(10));
1431
1432
1433   fcb->SetLineColor(1);
1434   fcb->SetNpx(1000);
1435   TLine* l1 = new TLine(m-n*s,0,m-n*s,fitTotal->GetParameter(4));
1436   TLine* l2 = new TLine(m+n*s,0,m+n*s,fitTotal->GetParameter(4));
1437   l1->SetLineColor(6);
1438   l2->SetLineColor(6);
1439   hfit->GetListOfFunctions()->Add(fcb);
1440   hfit->GetListOfFunctions()->Add(l1);
1441   hfit->GetListOfFunctions()->Add(l2);
1442
1443   Double_t cbParameters[7];
1444   Double_t covarianceMatrix[7][7];
1445
1446   for ( int ix = 0; ix < 7; ++ix )
1447   {
1448     cbParameters[ix] = fitTotal->GetParameter(ix+4);
1449   }
1450
1451   for ( int iy = 0; iy < 5; ++iy )
1452   {
1453     for ( int ix = 0; ix < 5; ++ix )
1454     {
1455       covarianceMatrix[ix][iy] = (fitResult->GetCovarianceMatrix())(ix+4,iy+4);
1456     }
1457   }
1458
1459   double njpsi = fcb->Integral(m-n*s,m+n*s)/h.GetBinWidth(1);
1460   double nerr = fcb->IntegralError(m-n*s,m+n*s,&cbParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(1);
1461
1462   r->Set("NofJpsi",njpsi,nerr);
1463
1464   return r;
1465 }
1466
1467 //_____________________________________________________________________________
1468 AliAnalysisMuMuResult* AliAnalysisMuMuResult::FitJpsi2CB2VWG(const TH1& h,
1469                                                              Double_t alphaLow,
1470                                                              Double_t nLow,
1471                                                              Double_t alphaUp,
1472                                                              Double_t nUp)
1473 {
1474   /// Fit using extended crystal ball + variable width gaussian
1475
1476   std::cout << Form("Fit with jpsi + psiprime VWG alphaLow=%5.2f nLow=%5.2f alphaUp=%5.2f nUp=%5.2f",
1477                     alphaLow,nLow,alphaUp,nUp) << std::endl;
1478
1479   Int_t nrebin = fMinv->GetXaxis()->GetNbins() / h.GetXaxis()->GetNbins();
1480
1481   TString resultName("JPSI2CB2VWG");
1482   if ( alphaLow > 0 )
1483   {
1484     resultName += TString::Format("alphaLow=%5.2f",alphaLow);
1485   }
1486   if ( nLow > 0 )
1487   {
1488     resultName += TString::Format("nLow=%5.2f",nLow);
1489   }
1490   if ( alphaUp > 0 )
1491   {
1492     resultName += TString::Format("alphaUp=%5.2f",alphaUp);
1493   }
1494   if ( nUp > 0 )
1495   {
1496     resultName += TString::Format("nUp=%5.2f",nUp);
1497   }
1498   resultName.ReplaceAll(" ","");
1499
1500   AliAnalysisMuMuResult* r = new AliAnalysisMuMuResult(h,resultName.Data(),nrebin);
1501
1502   TH1* hfit = r->Minv();
1503
1504   const Double_t xmin(2.2);
1505   const Double_t xmax(5.0);
1506
1507   TF1* fitTotal = new TF1("fitTotal",func2CB2VWG,xmin,xmax,12);
1508   fitTotal->SetParNames("kVWG","mVWG","sVWG1","sVWG2","kPsi","mPsi","sPsi","alPsi","nlPsi","auPsi","nuPsi");
1509   fitTotal->SetParName(11, "kPsi'");
1510
1511   fitTotal->SetParameter(0, 10000.);
1512   fitTotal->SetParameter(1, 1.9);
1513   fitTotal->SetParameter(2, 0.5);
1514   fitTotal->SetParLimits(2, 0., 100.);
1515   fitTotal->SetParameter(3, 0.3);
1516   fitTotal->SetParLimits(3, 0., 100.);
1517   fitTotal->SetParameter(4, 100.);
1518   fitTotal->SetParameter(5, 3.1);
1519   fitTotal->SetParLimits(5, 3.08, 3.2);
1520   fitTotal->SetParameter(6, 0.08);
1521   fitTotal->SetParLimits(6, 0.05, 0.15);
1522
1523 //  r = FitJpsi2CB2VWG(*hminv,0.93,5.59,2.32,3.39);
1524
1525   if ( alphaLow > 0 )
1526   {
1527     fitTotal->FixParameter(7, alphaLow);
1528   }
1529   else
1530   {
1531     fitTotal->SetParameter(7,0.9);
1532     fitTotal->SetParLimits(7,0.1,10.0);
1533   }
1534
1535   if ( nLow > 0 )
1536   {
1537     fitTotal->FixParameter(8, nLow);
1538   }
1539   else
1540   {
1541     fitTotal->SetParameter(8,5.0);
1542     fitTotal->SetParLimits(8,0.0,10.0);
1543   }
1544
1545   if ( alphaUp > 0 )
1546   {
1547     fitTotal->FixParameter(9, alphaUp);
1548   }
1549   else
1550   {
1551     fitTotal->SetParameter(9, 2.0);
1552     fitTotal->SetParLimits(9,0.1,10.0);
1553   }
1554
1555   if ( nUp > 0 )
1556   {
1557     fitTotal->FixParameter(10, nUp);
1558   }
1559   else
1560   {
1561     fitTotal->SetParameter(10,3.0);
1562     fitTotal->SetParLimits(10,0.0,10.0);
1563   }
1564
1565   fitTotal->SetParameter(11, 10.);
1566
1567   const char* fitOption = "SER"; //+";
1568
1569   TFitResultPtr fitResult = hfit->Fit(fitTotal,fitOption,"");
1570
1571   r->Set("MeanJpsi",fitTotal->GetParameter(5),fitTotal->GetParError(5));
1572   r->Set("SigmaJpsi",fitTotal->GetParameter(6),fitTotal->GetParError(6));
1573
1574   double m = r->GetValue("MeanJpsi");
1575   double s = r->GetValue("SigmaJpsi");
1576   double n = 3.0;
1577
1578   TF1* fcb = new TF1("fcb",CrystalBallExtended,xmin,xmax,7);
1579   fcb->SetParameters(fitTotal->GetParameter(4),
1580                      fitTotal->GetParameter(5),
1581                      fitTotal->GetParameter(6),
1582                      fitTotal->GetParameter(7),
1583                      fitTotal->GetParameter(8),
1584                      fitTotal->GetParameter(9),
1585                      fitTotal->GetParameter(10));
1586
1587
1588   fcb->SetLineColor(1);
1589   fcb->SetNpx(1000);
1590   TLine* l1 = new TLine(m-n*s,0,m-n*s,fitTotal->GetParameter(4));
1591   TLine* l2 = new TLine(m+n*s,0,m+n*s,fitTotal->GetParameter(4));
1592   l1->SetLineColor(6);
1593   l2->SetLineColor(6);
1594   hfit->GetListOfFunctions()->Add(fcb);
1595   hfit->GetListOfFunctions()->Add(l1);
1596   hfit->GetListOfFunctions()->Add(l2);
1597
1598   Double_t cbParameters[7];
1599   Double_t covarianceMatrix[7][7];
1600
1601   for ( int ix = 0; ix < 7; ++ix )
1602   {
1603     cbParameters[ix] = fitTotal->GetParameter(ix+4);
1604   }
1605
1606   for ( int iy = 0; iy < 5; ++iy )
1607   {
1608     for ( int ix = 0; ix < 5; ++ix )
1609     {
1610       covarianceMatrix[ix][iy] = (fitResult->GetCovarianceMatrix())(ix+4,iy+4);
1611     }
1612   }
1613
1614   double njpsi = fcb->Integral(m-n*s,m+n*s)/h.GetBinWidth(1);
1615   double nerr = fcb->IntegralError(m-n*s,m+n*s,&cbParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(1);
1616
1617   r->Set("NofJpsi",njpsi,nerr);
1618
1619   return r;
1620 }
1621
1622 //_____________________________________________________________________________
1623 AliAnalysisMuMuResult* AliAnalysisMuMuResult::FitJpsiNA48(const TH1& h)
1624 {
1625   /// fit using functional form from Ruben Shahoyan's thesis (2001) (eq. 4.8.)
1626
1627   std::cout << "Fit with jpsi NA50 Ruben eq. 4.8" << std::endl;
1628
1629   Int_t nrebin = fMinv->GetXaxis()->GetNbins() / h.GetXaxis()->GetNbins();
1630
1631   AliAnalysisMuMuResult* r = new AliAnalysisMuMuResult(h,"JPSINA",nrebin);
1632
1633   TH1* hfit = r->Minv();
1634
1635   const Double_t xmin(2.0);
1636   const Double_t xmax(5.0);
1637
1638   TF1* fitTotal = new TF1("fitTotal",funcJpsiNA48,xmin,xmax,12);
1639
1640   fitTotal->SetParName( 0, "c1");
1641   fitTotal->FixParameter(0,0.97);
1642
1643   fitTotal->SetParName( 1, "c2");
1644   fitTotal->FixParameter(1,1.05);
1645
1646   fitTotal->SetParName( 2, "d1");
1647   fitTotal->SetParameter(2,0.0);
1648   fitTotal->SetParLimits(2,0,1);
1649
1650   fitTotal->SetParName( 3, "d2");
1651   fitTotal->SetParameter(3,0.0);
1652   fitTotal->SetParLimits(3,0,1);
1653
1654   fitTotal->SetParName( 4, "g1");
1655   fitTotal->SetParameter(4,0.0);
1656   fitTotal->SetParLimits(4,0.0,1);
1657
1658   fitTotal->SetParName( 5, "g2");
1659   fitTotal->SetParameter(5,0.0);
1660   fitTotal->SetParLimits(5,0.0,1);
1661
1662   fitTotal->SetParName( 6, "m0");
1663   fitTotal->SetParameter(6,3.1);
1664   fitTotal->SetParLimits(6,2.8,3.4);
1665
1666   fitTotal->SetParName( 7, "sigma1");
1667   fitTotal->SetParameter(7,0.05);
1668   fitTotal->SetParLimits(7,0.01,0.2);
1669
1670   fitTotal->SetParName( 8, "sigma2");
1671   fitTotal->SetParameter(8,0.05);
1672   fitTotal->SetParLimits(8,0.01,0.2);
1673
1674   fitTotal->SetParName( 9, "b1");
1675   fitTotal->SetParameter(9,1.0);
1676   fitTotal->SetParLimits(9,0,1);
1677
1678   fitTotal->SetParName(10, "b2");
1679   fitTotal->SetParameter(10,1.0);
1680   fitTotal->SetParLimits(10,0,1);
1681
1682   fitTotal->SetParName(11, "norm");
1683   fitTotal->SetParameter(11,h.GetMaximum());
1684
1685   const char* fitOption = "SIER"; //+";
1686
1687   TFitResultPtr fitResult = hfit->Fit(fitTotal,fitOption,"");
1688
1689   r->Set("MeanJpsi",fitTotal->GetParameter(6),fitTotal->GetParError(6));
1690   r->Set("SigmaJpsi",
1691          fitTotal->GetParameter(7)+fitTotal->GetParameter(8),
1692          0.0);
1693
1694   double m = r->GetValue("MeanJpsi");
1695   double s = r->GetValue("SigmaJpsi");
1696   double n = 3.0;
1697
1698   TLine* l1 = new TLine(m-n*s,0,m-n*s,fitTotal->GetParameter(11));
1699   TLine* l2 = new TLine(m+n*s,0,m+n*s,fitTotal->GetParameter(11));
1700   l1->SetLineColor(6);
1701   l2->SetLineColor(6);
1702   hfit->GetListOfFunctions()->Add(l1);
1703   hfit->GetListOfFunctions()->Add(l2);
1704
1705   double njpsi = fitTotal->Integral(m-n*s,m+n*s)/h.GetBinWidth(1);
1706   double nerr = fitTotal->IntegralError(m-n*s,m+n*s)/h.GetBinWidth(1);
1707
1708   r->Set("NofJpsi",njpsi,nerr);
1709
1710   return r;
1711 }
1712
1713 /*
1714 //_____________________________________________________________________________
1715 void AliAnalysisMuMuResult::FitUpsilon(TH1& h)
1716 {
1717   std::cout << "Fit with upsilon alone" << std::endl;
1718
1719   const Double_t xmin(6.0);
1720   const Double_t xmax(12.0);
1721
1722   fitTotal = new TF1("fitTotal",funcCB2,xmin,xmax,7);
1723   fitTotal->SetParNames("N","alpha","n","mean","sigma","alphaprime","nprime");
1724   fitTotal->SetParameters(h.GetMaximum(),1,5,9.46,0.2,1.5,3);
1725   fitTotal->SetParLimits(0,0,h.GetMaximum()*2); // N
1726   fitTotal->SetParLimits(1,0,10); // alpha
1727   fitTotal->SetParLimits(2,1,10); // n
1728   fitTotal->SetParLimits(3,8,12); // mean
1729   fitTotal->SetParLimits(4,0.01,1); // sigma
1730   fitTotal->SetParLimits(5,0,10); // alpha
1731   fitTotal->SetParLimits(6,1,10); // n
1732
1733   h.Fit(fitTotal,"+","",6,12);
1734
1735
1736   Set("MeanUpsilon",fitTotal->GetParameter(3),fitTotal->GetParError(3));
1737   Set("SigmaUpsilon",fitTotal->GetParameter(4),fitTotal->GetParError(4));
1738
1739   double m = GetValue("MeanUpsilon");
1740   double s = GetValue("SigmaUpsilon");
1741   double n = 3.0;
1742
1743   Set("NofUpsilon",fitTotal->Integral(m-n*s,m+n*s)/h.GetBinWidth(1),fitTotal->IntegralError(m-n*s,m+n*s)/h.GetBinWidth(1));
1744 }
1745 */
1746
1747 //_____________________________________________________________________________
1748 Double_t AliAnalysisMuMuResult::ErrorAB(Double_t a, Double_t aerr, Double_t b, Double_t berr)
1749 {
1750   /// Compute the quadratic sum of 2 errors
1751
1752   Double_t e(0.0);
1753
1754   if ( TMath::Abs(a) > 1E-12 )
1755   {
1756     e += (aerr*aerr)/(a*a);
1757   }
1758
1759   if ( TMath::Abs(b) > 1E-12 )
1760   {
1761     e += (berr*berr)/(b*b);
1762   }
1763
1764   return TMath::Sqrt(e);
1765 }
1766
1767 //_____________________________________________________________________________
1768 Double_t AliAnalysisMuMuResult::ErrorABC(Double_t a, Double_t aerr, Double_t b, Double_t berr, Double_t c, Double_t cerror)
1769 {
1770   /// Compute the quadratic sum of 3 errors
1771
1772   Double_t e(0.0);
1773
1774   if ( TMath::Abs(a) > 1E-12 )
1775   {
1776     e += (aerr*aerr)/(a*a);
1777   }
1778
1779   if ( TMath::Abs(b) > 1E-12 )
1780   {
1781     e += (berr*berr)/(b*b);
1782   }
1783
1784   if ( TMath::Abs(b) > 1E-12 )
1785   {
1786     e += (cerror*cerror)/(c*c);
1787   }
1788
1789   return TMath::Sqrt(e);
1790 }
1791
1792
1793 //_____________________________________________________________________________
1794 Bool_t AliAnalysisMuMuResult::AddFit(const char* fitType, Int_t npar, Double_t* par)
1795 {
1796   // Add a fit to this result
1797
1798   TString msg(Form("fitType=%s npar=%d par[]=",fitType,npar));
1799
1800   for ( Int_t i = 0; i < npar; ++i )
1801   {
1802     msg += TString::Format("%e,",par[i]);
1803   }
1804
1805   msg += TString::Format(" minv=%p %d",fMinv,fMinv?TMath::Nint(fMinv->GetEntries()):0);
1806
1807   if ( !fMinv ) return kFALSE;
1808
1809   TString sFitType(fitType);
1810   sFitType.ToUpper();
1811   Int_t nrebin(1);
1812
1813   if (sFitType.CountChar(':'))
1814   {
1815     Int_t index = sFitType.Index(":");
1816     nrebin = TString(sFitType(index+1,sFitType.Length()-index-1)).Atoi();
1817     sFitType = sFitType(0,index);
1818   }
1819
1820   msg += TString::Format(" nrebin=%d",nrebin);
1821
1822   AliDebug(1,msg.Data());
1823
1824
1825   if ( fMinv->GetEntries()<100 && !sFitType.Contains("COUNT")) return kFALSE;
1826
1827   TH1* hminv = static_cast<TH1*>(fMinv->Clone());
1828
1829   hminv->Rebin(nrebin);
1830   hminv->SetDirectory(0);
1831
1832   AliAnalysisMuMuResult* r(0x0);
1833
1834   if ( sFitType=="PSI1")
1835   {
1836     r = FitJpsi(*hminv);
1837   }
1838   else if ( sFitType == "PSILOW")
1839   {
1840     r = FitJpsi2CB2VWG(*hminv,-1,-1,-1,-1); // free tails
1841   }
1842   else if ( sFitType == "PSILOWMCTAILS" )
1843   {
1844     if ( npar!= 4 )
1845     {
1846       AliError("Cannot use PSILOWMCTAILS without being given the MC tails !");
1847       delete hminv;
1848       return kFALSE;
1849     }
1850     r = FitJpsi2CB2VWG(*hminv,par[0],par[1],par[2],par[3]);
1851     if (r)
1852     {
1853       r->SetAlias("MCTAILS");
1854     }
1855   }
1856   else if ( sFitType.BeginsWith("PSILOWALPHA") )
1857   {
1858     Float_t lpar[] = { 0.0, 0.0, 0.0, 0.0 };
1859
1860     AliDebug(1,Form("sFitType=%s",sFitType.Data()));
1861
1862     sscanf(sFitType.Data(),"PSILOWALPHALOW%fNLOW%fALPHAUP%fNUP%f",
1863            &lpar[0],&lpar[1],&lpar[2],&lpar[3]);
1864
1865     AliDebug(1,Form("PSILOW ALPHALOW=%f NLOW=%f ALPHAUP=%f NUP=%f",lpar[0],lpar[1],lpar[2],lpar[3]));
1866
1867     if ( lpar[0] == 0.0 && lpar[1] == 0.0 && lpar[0] == 0.0 && lpar[1] == 0.0 )
1868     {
1869       AliError("Cannot work with zero tails !");
1870     }
1871     else
1872     {
1873       r = FitJpsi2CB2VWG(*hminv,lpar[0],lpar[1],lpar[2],lpar[3]);
1874     }
1875   }
1876   else if ( sFitType == "COUNTJPSI" )
1877   {
1878     r = new AliAnalysisMuMuResult(*hminv,"COUNTJPSI",1);
1879     Double_t n = CountParticle(*hminv,"Jpsi");
1880     r->Set("NofJpsi",n,TMath::Sqrt(n));
1881   }
1882
1883
1884   if ( r )
1885   {
1886     r->Print();
1887     r->SetBin(Bin());
1888     r->SetNofTriggers(NofTriggers());
1889     r->SetNofRuns(NofRuns());
1890
1891     if (!fSubResults)
1892     {
1893       fSubResults = new TObjArray;
1894       fSubResults->SetOwner(kTRUE);
1895     }
1896
1897     fSubResults->Add(r);
1898   }
1899
1900   delete hminv;
1901
1902   return (r!=0x0);
1903 }
1904
1905 //_____________________________________________________________________________
1906 Double_t AliAnalysisMuMuResult::GetErrorStat(const char* name, const char* subResultName) const
1907 {
1908   // compute the mean value from all subresults
1909
1910   if ( strlen(subResultName) > 0 )
1911   {
1912     if ( !fSubResults)
1913     {
1914       AliError(Form("No subresult from which I could get the %s one...",subResultName));
1915       return TMath::Limits<Double_t>::Max();
1916     }
1917     AliAnalysisMuMuResult* sub = static_cast<AliAnalysisMuMuResult*>(fSubResults->FindObject(subResultName));
1918     if (!sub)
1919     {
1920       AliError(Form("Could not get subresult named %s",subResultName));
1921       return TMath::Limits<Double_t>::Max();
1922     }
1923     return sub->GetErrorStat(name);
1924   }
1925
1926   if ( fMap )
1927   {
1928     TObjArray* p = static_cast<TObjArray*>(fMap->GetValue(name));
1929     if (p)
1930     {
1931       TParameter<double>* val = static_cast<TParameter<double>*>(p->At(kErrorStat));
1932       return val->GetVal();
1933     }
1934   }
1935
1936   TIter next(fSubResults);
1937   AliAnalysisMuMuResult* r;
1938   Int_t n(0);
1939   Double_t mean(0);
1940
1941   while ( ( r = static_cast<AliAnalysisMuMuResult*>(next()) ) )
1942   {
1943     if ( r->HasValue(name) )
1944     {
1945       mean += r->GetErrorStat(name);
1946       ++n;
1947     }
1948   }
1949   return ( n ? mean/n : 0.0 );
1950 }
1951
1952 //_____________________________________________________________________________
1953 Double_t AliAnalysisMuMuResult::GetValue(const char* name, const char* subResultName) const
1954 {
1955   // get a value (either directly or by computing the mean of the subresults)
1956
1957   if ( strlen(subResultName) > 0 )
1958   {
1959     if ( !fSubResults)
1960     {
1961       AliError(Form("No subresult from which I could get the %s one...",subResultName));
1962       return TMath::Limits<Double_t>::Max();
1963     }
1964     AliAnalysisMuMuResult* sub = static_cast<AliAnalysisMuMuResult*>(fSubResults->FindObject(subResultName));
1965     if (!sub)
1966     {
1967       AliError(Form("Could not get subresult named %s",subResultName));
1968       return TMath::Limits<Double_t>::Max();
1969     }
1970     return sub->GetValue(name);
1971   }
1972
1973   if (fMap)
1974   {
1975     TObjArray* p = static_cast<TObjArray*>(fMap->GetValue(name));
1976     if (p)
1977     {
1978       TParameter<double>* val = static_cast<TParameter<double>*>(p->At(kValue));
1979       return val->GetVal();
1980     }
1981   }
1982
1983   // compute the mean value from all subresults
1984   TIter next(fSubResults);
1985   AliAnalysisMuMuResult* r;
1986   Int_t n(0);
1987   Double_t mean(0);
1988
1989   while ( ( r = static_cast<AliAnalysisMuMuResult*>(next()) ) )
1990   {
1991     if ( r->HasValue(name) )
1992     {
1993       mean += r->GetValue(name);
1994       ++n;
1995     }
1996   }
1997   return ( n ? mean/n : 0.0 );
1998 }
1999
2000 //_____________________________________________________________________________
2001 Bool_t AliAnalysisMuMuResult::HasValue(const char* name, const char* subResultName) const
2002 {
2003   /// Whether this result (or subresult if subResultName is provided) has a property
2004   /// named "name"
2005
2006   if ( strlen(subResultName) > 0 )
2007   {
2008     if ( !fSubResults)
2009     {
2010       AliError(Form("No subresult from which I could get the %s one...",subResultName));
2011       return kFALSE;
2012     }
2013     AliAnalysisMuMuResult* sub = static_cast<AliAnalysisMuMuResult*>(fSubResults->FindObject(subResultName));
2014     if (!sub)
2015     {
2016       AliError(Form("Could not get subresult named %s",subResultName));
2017       return kFALSE;
2018     }
2019     return sub->HasValue(name);
2020   }
2021
2022   if ( fMap && ( fMap->GetValue(name) != 0x0 ) )
2023   {
2024     return kTRUE;
2025   }
2026
2027   TIter next(fSubResults);
2028   AliAnalysisMuMuResult* r;
2029
2030   while ( ( r = static_cast<AliAnalysisMuMuResult*>(next()) ) )
2031   {
2032     if ( r->HasValue(name) ) return kTRUE;
2033   }
2034
2035   return kFALSE;
2036 }
2037
2038 //_____________________________________________________________________________
2039 THashList* AliAnalysisMuMuResult::Keys() const
2040 {
2041   /// Return the complete list of keys we're using
2042   if (!fKeys)
2043   {
2044     fKeys = new THashList;
2045     fKeys->SetOwner(kTRUE);
2046     TIter next(fMap);
2047     TObjString* key;
2048
2049     while ( ( key = static_cast<TObjString*>(next()) ) )
2050     {
2051       if ( !fKeys->FindObject(key->String()) )
2052       {
2053         fKeys->Add(new TObjString(key->String()));
2054       }
2055     }
2056
2057     AliAnalysisMuMuResult* r;
2058     TIter nextResult(fSubResults);
2059
2060     while ( ( r = static_cast<AliAnalysisMuMuResult*>(nextResult())) )
2061     {
2062       TIter nextHL(r->Keys());
2063       TObjString* s;
2064
2065       while ( ( s = static_cast<TObjString*>(nextHL())) )
2066       {
2067         if ( !fKeys->FindObject(s->String()) )
2068         {
2069           fKeys->Add(new TObjString(s->String()));
2070         }
2071       }
2072     }
2073
2074   }
2075   return fKeys;
2076 }
2077
2078 //_____________________________________________________________________________
2079 Long64_t AliAnalysisMuMuResult::Merge(TCollection* list)
2080 {
2081   /// Merge method
2082   ///
2083   /// Merge a list of AliAnalysisMuMuResult objects with this
2084   /// Returns the number of merged objects (including this).
2085   ///
2086   /// Note that the merging is to be understood here as an average operation
2087
2088   AliInfo(Form("this=%p",this));
2089   if (!list) return 0;
2090
2091   if (list->IsEmpty()) return 1;
2092
2093   TIter next(list);
2094   TObject* currObj;
2095   TList hminvList;
2096   hminvList.SetOwner(kTRUE);
2097
2098   Double_t thisWeight = Weight();
2099   Double_t sumOfWeights = thisWeight;
2100
2101   Double_t nofTriggers = fNofTriggers*thisWeight;
2102   Double_t nofRuns = fNofRuns*thisWeight;
2103
2104   fNofRuns *= thisWeight;
2105
2106   while ( ( currObj = next() ) )
2107   {
2108     AliAnalysisMuMuResult* result = dynamic_cast<AliAnalysisMuMuResult*>(currObj);
2109     if (!result)
2110     {
2111       AliFatal(Form("object named \"%s\" is a %s instead of an AliAnalysisMuMuResult!", currObj->GetName(), currObj->ClassName()));
2112       continue;
2113     }
2114
2115     Double_t w = result->Weight();
2116
2117     nofRuns += result->NofRuns()*w;
2118     nofTriggers += result->NofTriggers()*w;
2119     fWeight += result->fWeight;
2120     sumOfWeights += w;
2121   }
2122
2123   thisWeight/= sumOfWeights;
2124   fNofRuns = nofRuns/sumOfWeights;
2125   fNofTriggers = nofTriggers/sumOfWeights;
2126   fWeight /= sumOfWeights;
2127
2128   AliInfo(Form("thisWeight=%e sumOfWeight=%8.2f noftriggers=%e weight=%e",thisWeight,sumOfWeights,1.0*fNofTriggers,fWeight));
2129
2130   TIter nextKey(fMap);
2131   TObjString* key;
2132
2133   while ( ( key = static_cast<TObjString*>(nextKey())) )
2134   {
2135     AliInfo(key->String().Data());
2136
2137     Double_t value = GetValue(key->String())*thisWeight;
2138     Double_t estat = GetErrorStat(key->String())*GetErrorStat(key->String())*thisWeight*thisWeight;
2139
2140     Double_t test(thisWeight);
2141
2142     next.Reset();
2143
2144     while ( ( currObj = next() ) )
2145     {
2146       AliAnalysisMuMuResult* result = dynamic_cast<AliAnalysisMuMuResult*>(currObj);
2147
2148       if (!result)
2149       {
2150         continue;
2151       }
2152
2153       if (!result->HasValue(key->String()))
2154       {
2155         AliError(Form("Did not find key %s in of the result to merge",key->String().Data()));
2156         continue;
2157       }
2158
2159       // can only merge under the condition we have the same bin
2160
2161       if ( fBin != result->Bin() )
2162       {
2163         AliError("Cannot merge results of different bin");
2164         continue;
2165       }
2166
2167       Double_t w = result->Weight()/sumOfWeights;
2168
2169       Double_t w2 = w*w;
2170
2171       test += w;
2172
2173       value += result->GetValue(key->String())*w;
2174       estat += result->GetErrorStat(key->String())*result->GetErrorStat(key->String())*w2;
2175
2176     }
2177
2178     Set(key->String(),
2179         value,
2180         TMath::Sqrt(estat));
2181
2182     AliInfo(Form("test=%e",test));
2183   }
2184
2185   if ( fMinv )
2186   {
2187     fMinv->Scale(thisWeight);
2188   }
2189
2190   next.Reset();
2191
2192   while ( ( currObj = next() ) )
2193   {
2194     AliAnalysisMuMuResult* result = dynamic_cast<AliAnalysisMuMuResult*>(currObj);
2195
2196     if ( result->Minv() )
2197     {
2198       TH1* h = static_cast<TH1*>(result->Minv()->Clone());
2199       AliInfo(Form("Nbins %d xmin %e xmax %e",h->GetXaxis()->GetNbins(),h->GetXaxis()->GetXmin(),
2200                    h->GetXaxis()->GetXmax()));
2201       h->Scale(result->Weight());
2202       hminvList.Add(h);
2203     }
2204   }
2205
2206   if ( fMinv )
2207   {
2208     fMinv->Merge(&hminvList);
2209     fMinv->Scale(1.0/sumOfWeights);
2210   }
2211
2212   TIter nextSubresult(fSubResults);
2213   AliAnalysisMuMuResult* r;
2214
2215   while ( ( r = static_cast<AliAnalysisMuMuResult*>(nextSubresult())) )
2216   {
2217     TList sublist;
2218
2219     next.Reset();
2220
2221     while ( ( currObj = next() ) )
2222     {
2223       sublist.Add(currObj);
2224     }
2225
2226     r->Merge(&sublist);
2227   }
2228
2229   return list->GetEntries()+1;
2230 }
2231
2232 //_____________________________________________________________________________
2233 Int_t AliAnalysisMuMuResult::NofRuns() const
2234 {
2235   /// Get the number of runs
2236   if ( !Mother() ) return fNofRuns;
2237   else return Mother()->NofRuns();
2238 }
2239
2240 //_____________________________________________________________________________
2241 Int_t AliAnalysisMuMuResult::NofTriggers() const
2242 {
2243   /// Get the number of triggers
2244
2245   if ( !Mother() ) return fNofTriggers;
2246   else return Mother()->NofTriggers();
2247 }
2248
2249 //_____________________________________________________________________________
2250 void AliAnalysisMuMuResult::Print(Option_t* opt) const
2251 {
2252   /// printout
2253
2254   TString sopt(opt);
2255   sopt.ToUpper();
2256
2257   for ( Int_t i = 0; i < 9; ++i )
2258   {
2259     sopt.ReplaceAll(Form("%d",i),"");
2260   }
2261
2262   TString pot(sopt);
2263   pot.ReplaceAll("ALL","");
2264   pot.ReplaceAll("FULL","");
2265
2266   std::cout << pot.Data();
2267
2268   if ( fAlias.Length() > 0 )
2269   {
2270     std::cout << Form("%s - ",fAlias.Data());
2271   }
2272
2273   std::cout << Form("%50s - NRUNS %d - NTRIGGER %10d - %s%s",
2274                GetName(),
2275                     NofRuns(),
2276                     NofTriggers(),
2277                                   fWeight > 0.0 ? Form(" WEIGHT %e -",fWeight) : "",
2278                            fBin.AsString().Data());
2279
2280     if ( fSubResults && fSubResults->GetEntries()>1 )
2281     {
2282       std::cout << " (" << fSubResults->GetEntries()-1 << " subresults)";
2283     }
2284
2285   if (!fSubResults )
2286   {
2287     std::cout << Form(" - REBIN %d",fRebin) << std::endl;
2288   }
2289
2290   std::cout << std::endl;
2291
2292   if ( sopt.Contains("DUMP") )
2293   {
2294     TIter next(fMap);
2295     TObjString* str;
2296     while ( ( str = static_cast<TObjString*>(next()) ) )
2297     {
2298       TObjArray* a = static_cast<TObjArray*>(fMap->GetValue(str->String()));
2299
2300       std::cout << Form("%s %e %e",
2301                         str->String().Data(),
2302                         static_cast<TParameter<Double_t>*> (a->At(kValue))->GetVal(),
2303                         static_cast<TParameter<Double_t>*> (a->At(kErrorStat))->GetVal()) << std::endl;
2304     }
2305   }
2306   else
2307   {
2308
2309     PrintParticle("Jpsi",pot.Data());
2310     PrintParticle("PsiPrime",pot.Data());
2311     PrintParticle("Upsilon",pot.Data());
2312   }
2313
2314   if ( HasValue("MBR"))
2315   {
2316     std::cout << Form("\t\tMBR %e +- %e",GetValue("MBR"),GetErrorStat("MBR")) << std::endl;
2317   }
2318
2319   if ( fSubResults /* && fSubResults->GetEntries() > 1 */ && ( sopt.Contains("ALL") || sopt.Contains("FULL") ) )
2320   {
2321     std::cout << pot.Data() << "\t===== sub results ===== " << std::endl;
2322
2323     sopt += "\t\t";
2324
2325     TIter next(fSubResults);
2326     AliAnalysisMuMuResult* r;
2327
2328     while ( ( r = static_cast<AliAnalysisMuMuResult*>(next()) ) )
2329     {
2330       r->Print(sopt.Data());
2331     }
2332   }
2333 }
2334
2335 //_____________________________________________________________________________
2336 void AliAnalysisMuMuResult::PrintParticle(const char* particle, const char* opt) const
2337 {
2338   /// Print all information about one particule type
2339
2340   Double_t npart = GetValue(Form("Nof%s",particle));
2341   if (npart<=0) return;
2342
2343
2344   std::cout << opt << Form("\t%s",particle) << std::endl;
2345
2346   //  Double_t npartError = GetErrorStat(Form("Nof%s",particle));
2347 //  std::cout << opt << Form("\t\t%20s %9.2f +- %5.2f","Count",npart,npartError) << std::endl;
2348
2349   TIter next(Keys());
2350   TObjString* key;
2351
2352   while ( ( key = static_cast<TObjString*>(next()) ) )
2353   {
2354     if ( !key->String().Contains(particle) ) continue;
2355
2356     PrintValue(key->String(),opt,GetValue(key->String()),GetErrorStat(key->String()));
2357   }
2358 }
2359
2360 //_____________________________________________________________________________
2361 void AliAnalysisMuMuResult::PrintValue(const char* key, const char* opt, Double_t value, Double_t errorStat)
2362 {
2363   // print one value and its associated error
2364
2365   if ( TString(key).Contains("AccEff") )
2366   {
2367     std::cout << opt << Form("\t\t%20s %9.2f +- %5.2f %%",key,value*100,errorStat*100) << std::endl;
2368   }
2369   else if ( TString(key).BeginsWith("Sigma") ||  TString(key).BeginsWith("Mean") )
2370   {
2371     std::cout << opt << Form("\t\t%20s %9.2f +- %5.2f MeV/c^2",key,value*1E3,1E3*errorStat) << std::endl;
2372   }
2373   else if ( TString(key).Contains("Nof") )
2374   {
2375     std::cout << opt << Form("\t\t%20s %9.2f +- %5.2f",key,value,errorStat) << std::endl;
2376   }
2377   else if ( value > 1E-3 && value < 1E3 )
2378   {
2379     std::cout << opt << Form("\t\t%20s %9.2f +- %5.2f",key,value,errorStat) << std::endl;
2380   }
2381   else
2382   {
2383     std::cout << opt << Form("\t\t%20s %9.2e +- %9.2e",key,value,errorStat) << std::endl;
2384   }
2385 }
2386
2387 //_____________________________________________________________________________
2388 void AliAnalysisMuMuResult::Set(const char* name, Double_t value, Double_t errorStat)
2389 {
2390   /// Set a (value,error) pair with a given name
2391
2392   if ( !fMap )
2393   {
2394     fMap = new TMap;
2395     fMap->SetOwnerKeyValue(kTRUE,kTRUE);
2396   }
2397
2398   TObjArray* p = static_cast<TObjArray*>(fMap->GetValue(name));
2399   if (!p)
2400   {
2401     p = new TObjArray(4);
2402
2403     p->SetOwner(kTRUE);
2404
2405     p->AddAt(new TParameter<Double_t>(name,value),kValue);
2406     p->AddAt(new TParameter<Double_t>(name,errorStat),kErrorStat);
2407
2408     fMap->Add(new TObjString(name),p);
2409   }
2410   else
2411   {
2412     static_cast<TParameter<double>*>(p->At(kValue))->SetVal(value);
2413     static_cast<TParameter<double>*>(p->At(kErrorStat))->SetVal(errorStat);
2414   }
2415
2416   if ( TString(name)=="NofJpsi" )
2417   {
2418     if ( NofTriggers() > 0 )
2419     {
2420       Double_t rate = value/NofTriggers();
2421       Double_t rateError = rate*ErrorAB(value,errorStat,NofTriggers(),TMath::Sqrt(NofTriggers()));
2422       Set("RateJpsi",rate,rateError);
2423     }
2424   }
2425 }
2426
2427 //_____________________________________________________________________________
2428 void AliAnalysisMuMuResult::SetBin(const AliAnalysisMuMuBinning::Range& bin)
2429 {
2430   /// Set the bin
2431
2432   if (!Mother()) fBin = bin;
2433   else Mother()->SetBin(bin);
2434 }
2435
2436 //_____________________________________________________________________________
2437 void AliAnalysisMuMuResult::SetNofInputParticles(const TH1& hminv)
2438 {
2439   /// Set the number of input particle from the invariant mass spectra
2440
2441   static const char* particleNames[] = { "Jpsi" , "PsiPrime", "Upsilon","UpsilonPrime" };
2442
2443   for ( Int_t i = 0; i < 4; ++i )
2444   {
2445     Double_t n = CountParticle(hminv,particleNames[i]);
2446
2447     AliDebug(1,Form("i=%d particle %s n %e",i,particleNames[i],n));
2448
2449     if ( n > 0 )
2450     {
2451       SetNofInputParticles(particleNames[i],n);
2452     }
2453   }
2454 }
2455
2456 //_____________________________________________________________________________
2457 void AliAnalysisMuMuResult::SetNofInputParticles(const char* particle, int n)
2458 {
2459   /// Set the number of input particles (so it is a MC result)
2460   /// and (re)compute the AccxEff values
2461
2462   Set(Form("NofInput%s",particle),n,TMath::Sqrt(n));
2463
2464   if (n<=0)
2465   {
2466     Set(Form("AccEff%s",particle),0,0);
2467     return;
2468   }
2469
2470   Double_t npart = GetValue(Form("Nof%s",particle));
2471   Double_t npartErr  = GetErrorStat(Form("Nof%s",particle));
2472   Double_t ninput = GetValue(Form("NofInput%s",particle));
2473   Double_t ninputErr = GetErrorStat(Form("NofInput%s",particle));
2474
2475   Set(Form("AccEff%s",particle),
2476       npart/ninput,
2477       (npart/ninput)*ErrorAB(npart,npartErr,ninput,ninputErr));
2478
2479   TIter next(fSubResults);
2480   AliAnalysisMuMuResult* r;
2481
2482   while ( ( r = static_cast<AliAnalysisMuMuResult*>(next())) )
2483   {
2484     r->Set(Form("NofInput%s",particle),n,TMath::Sqrt(n));
2485
2486     npart = r->GetValue(Form("Nof%s",particle));
2487     npartErr = r->GetErrorStat(Form("Nof%s",particle));
2488
2489     r->Set(Form("AccEff%s",particle),
2490            npart/ninput,
2491            (npart/ninput)*ErrorAB(npart,npartErr,ninput,ninputErr));
2492
2493   }
2494 }
2495
2496 //_____________________________________________________________________________
2497 void AliAnalysisMuMuResult::SetNofRuns(Int_t n)
2498 {
2499   if ( !Mother() ) fNofRuns=n;
2500   else Mother()->SetNofRuns(n);
2501 }
2502
2503 //_____________________________________________________________________________
2504 void AliAnalysisMuMuResult::SetNofTriggers(Int_t n)
2505 {
2506   if ( !Mother() ) fNofTriggers=n;
2507   else Mother()->SetNofTriggers(n);
2508 }
2509
2510 //_____________________________________________________________________________
2511 void AliAnalysisMuMuResult::SetMinv(const TH1& hminv)
2512 {
2513     /// Set the inv. mass spectrum to be fitted.
2514   static UInt_t n(0);
2515
2516   delete fMinv;
2517   fMinv = static_cast<TH1*>(hminv.Clone(Form("Minv%u",n++)));
2518   fMinv->SetDirectory(0);
2519 }
2520
2521 //_____________________________________________________________________________
2522 AliAnalysisMuMuResult*
2523 AliAnalysisMuMuResult::SubResult(const char* subResultName) const
2524 {
2525   /// get a given subresult
2526   if (!fSubResults)
2527   {
2528     return 0x0;
2529   }
2530   TIter next(fSubResults);
2531   AliAnalysisMuMuResult* r;
2532   while ( ( r = static_cast<AliAnalysisMuMuResult*>(next())) )
2533   {
2534     if ( r->Alias() == subResultName )
2535     {
2536       return r;
2537     }
2538   }
2539   return 0x0;
2540 }
2541