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