PWG/muondep/AliAnalysisMuMuJpsiResult.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 /// Class to hold results about J/psi
  18 /// like number of of J/psi (before and after Acc x Eff correction),
  19 /// Acc x Eff correction, Yield, RAB, etc...
  20 ///
  21 /// author: Laurent Aphecetche (Subatech)
  22 ///
  23
  24 #include "AliAnalysisMuMuJpsiResult.h"
  25
  26 ClassImp(AliAnalysisMuMuJpsiResult)
  27
  28 #include "TF1.h"
  29 #include "TFitResult.h"
  30 #include "TH1.h"
  31 #include "TH2.h"
  32 #include "THashList.h"
  33 #include "TLine.h"
  34 #include "TList.h"
  35 #include "TMap.h"
  36 #include "TMath.h"
  37 #include "TObjArray.h"
  38 #include "TParameter.h"
  39 #include "AliAnalysisMuMuBinning.h"
  40 #include "AliLog.h"
  41 #include <map>
  42
  43 namespace {
  44
  45   const std::map<std::string,Double_t>& MassMap()
  46   {
  47     /// a simple map of masses...
  48     static std::map<std::string,Double_t> massMap;
  49     // not super elegant, but that way we do not depend on TDatabasePDG and thus
  50     // can decide on our particle naming
  51     if (massMap.empty())
  52     {
  53       massMap["Jpsi"] = 3.096916e+00;
  54       massMap["PsiPrime"] = 3.68609e+00;
  55       massMap["Upsilon"] = 9.46030e+00;
  56       massMap["UpsilonPrime"] = 1.00233e+01;
  57     }
  58     return massMap;
  59   }
  60
  61
  62   Double_t funcCB(Double_t* xx, Double_t* par)
  63   {
  64     /// Crystal ball
  65
  66     Double_t norm = par[0];
  67     Double_t alpha = par[1];
  68     Double_t n = par[2];
  69     Double_t mean = par[3];
  70     Double_t sigma = par[4];
  71
  72     Double_t x = xx[0];
  73
  74     Double_t a = TMath::Power(n/TMath::Abs(alpha),n)*TMath::Exp(-0.5*alpha*alpha);
  75     Double_t b = n/TMath::Abs(alpha) - TMath::Abs(alpha);
  76
  77     Double_t y = ( TMath::Abs(sigma) > 1E-12 ? (x-mean)/sigma : 0 );
  78
  79     if ( y > alpha*-1.0 )
  80     {
  81       return norm*TMath::Exp(-0.5*y*y);
  82     }
  83     else
  84     {
  85       return norm*a*TMath::Power(b-y,-n);
  86     }
  87   }
  88
  89   Double_t funcJpsiGCBE(Double_t* xx, Double_t* par)
  90   {
  91     /// crystal ball + expo + gaussian
  92     Double_t x = xx[0];
  93
  94     Double_t g = par[0]*TMath::Gaus(x,par[1],par[2]);
  95
  96     Double_t jpsi = funcCB(xx,par+3);
  97
  98     Double_t expo = par[8]*TMath::Exp(par[9]*x);
  99
 100     return g+expo+jpsi;
 101   }
 102
 103   Double_t funcCB2(Double_t* xx, Double_t* par)
 104   {
 105     /// CB2 = extended crystal ball
 106
 107     Double_t norm = par[0];
 108     Double_t alpha = par[1];
 109     Double_t n = par[2];
 110     Double_t mean = par[3];
 111     Double_t sigma = par[4];
 112     Double_t alphaprime = par[5];
 113     Double_t nprime = par[6];
 114
 115     Double_t x = xx[0];
 116
 117     Double_t a = TMath::Power(n/TMath::Abs(alpha),n)*TMath::Exp(-0.5*alpha*alpha);
 118     Double_t b = n/TMath::Abs(alpha) - TMath::Abs(alpha);
 119     Double_t c = TMath::Power(nprime/TMath::Abs(alphaprime),nprime)*TMath::Exp(-0.5*alphaprime*alphaprime);
 120     Double_t d = nprime/TMath::Abs(alphaprime) - TMath::Abs(alphaprime);
 121
 122     Double_t y = ( TMath::Abs(sigma) > 1E-12 ? (x-mean)/sigma : 0 );
 123
 124     if ( y > alphaprime )
 125     {
 126       return norm*c*TMath::Power(d+y,-nprime);
 127     }
 128     else if ( y > alpha*-1.0 )
 129     {
 130       return norm*TMath::Exp(-0.5*y*y);
 131     }
 132     else
 133     {
 134       return norm*a*TMath::Power(b-y,-n);
 135     }
 136   }
 137
 138
 139   Double_t funcJpsiNA48(Double_t*x, Double_t* par)
 140   {
 141     /// Fit function from e.q. 4.8 of Ruben's PhD.
 142     Double_t c1 = par[0];
 143     Double_t c2 = par[1];
 144     Double_t d1 = par[2];
 145     Double_t d2 = par[3];
 146     Double_t g1 = par[4];
 147     Double_t g2 = par[5];
 148     Double_t m0 = par[6];
 149     Double_t sigma1 = par[7];
 150     Double_t sigma2 = par[8];
 151     Double_t b1 = par[9];
 152     Double_t b2 = par[10];
 153     Double_t norm = par[11];
 154
 155     Double_t m = x[0];
 156
 157     Double_t rv(0);
 158
 159     if ( m <= c1*m0 )
 160     {
 161       Double_t e = d1-g1*TMath::Sqrt(c1*m0-m);
 162       rv = TMath::Power(b1*(c1*m0-m),e);
 163       rv += sigma1;
 164     }
 165     else if( m >= c1*m0 && m <= m0 )
 166     {
 167       rv = sigma1;
 168     }
 169     else if ( m >= m0 && m < c2*m0 )
 170     {
 171       rv = sigma2;
 172     }
 173     else if( m >= c2*m0 )
 174     {
 175       Double_t e = d2-g2*TMath::Sqrt(m-c2*m0);
 176       rv = TMath::Power(b2*(m-c2*m0),e);
 177       rv += sigma2;
 178     }
 179
 180     return norm*TMath::Exp(-(m-m0)*(m-m0)/(2.0*rv*rv));
 181   }
 182
 183   //------------------------------------------------------------------------------
 184   Double_t BackgroundVWG(Double_t *x, Double_t *par)
 185   {
 186     // gaussian variable width
 187     Double_t sigma = par[2]+par[3]*((x[0]-par[1])/par[1]);
 188     return par[0]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/(2.*sigma*sigma));
 189
 190   }
 191
 192   //------------------------------------------------------------------------------
 193   Double_t CrystalBallExtended(Double_t *x,Double_t *par)
 194   {
 195     //par[0] = Normalization
 196     //par[1] = mean
 197     //par[2] = sigma
 198     //par[3] = alpha
 199     //par[4] = n
 200     //par[5] = alpha'
 201     //par[6] = n'
 202
 203     Double_t t = (x[0]-par[1])/par[2];
 204     if (par[3] < 0) t = -t;
 205
 206     Double_t absAlpha = fabs((Double_t)par[3]);
 207     Double_t absAlpha2 = fabs((Double_t)par[5]);
 208
 209     if (t >= -absAlpha && t < absAlpha2) // gaussian core
 210     {
 211       return par[0]*(exp(-0.5*t*t));
 212     }
 213
 214     if (t < -absAlpha) //left tail
 215     {
 216       Double_t a =  TMath::Power(par[4]/absAlpha,par[4])*exp(-0.5*absAlpha*absAlpha);
 217       Double_t b = par[4]/absAlpha - absAlpha;
 218       return par[0]*(a/TMath::Power(b - t, par[4]));
 219     }
 220
 221     if (t >= absAlpha2) //right tail
 222     {
 223
 224       Double_t c =  TMath::Power(par[6]/absAlpha2,par[6])*exp(-0.5*absAlpha2*absAlpha2);
 225       Double_t d = par[6]/absAlpha2 - absAlpha2;
 226       return par[0]*(c/TMath::Power(d + t, par[6]));
 227     }
 228
 229     return 0. ;
 230   }
 231
 232
 233   //---------------------------------------------------------------------------
 234 //  Double_t fitFunctionVWG(Double_t *x, Double_t *par)
 235 //  {
 236 //    if (x[0] > 2.9 && x[0] < 3.3) TF1::RejectPoint();
 237 //    return BackgroundVWG(x, par);
 238 //  }
 239
 240   //---------------------------------------------------------------------------
 241   Double_t fitFunctionCB2VWG(Double_t *x, Double_t *par)
 242   {
 243     return BackgroundVWG(x, par) + CrystalBallExtended(x, &par[4]);
 244   }
 245
 246   //---------------------------------------------------------------------------
 247   Double_t func2CB2VWG(Double_t *x, Double_t *par)
 248   {
 249     /// 2 extended crystal balls + variable width gaussian
 250     /// width of the second CB related to the first (free) one.
 251
 252     Double_t par2[7] = {
 253       par[11],
 254       3.68609+(par[5]-3.096916)/3.096916*3.68609,
 255       par[6]/3.096916*3.68609,
 256       par[7],
 257       par[8],
 258       par[9],
 259       par[10]
 260     };
 261     return BackgroundVWG(x, par) + CrystalBallExtended(x, &par[4]) + CrystalBallExtended(x, par2);
 262   }
 263 }
 264
 265 //_____________________________________________________________________________
 266 AliAnalysisMuMuJpsiResult::AliAnalysisMuMuJpsiResult(TRootIOCtor* /*io*/) :
 267 AliAnalysisMuMuResult("",""),
 268 fNofRuns(),
 269 fNofTriggers(-1),
 270 fMinv(0x0),
 271 fBin(),
 272 fRebin(0),
 273 fTriggerClass(),
 274 fEventSelection(),
 275 fPairSelection(),
 276 fCentralitySelection()
 277 {
 278 }
 279
 280 //_____________________________________________________________________________
 281 AliAnalysisMuMuJpsiResult::AliAnalysisMuMuJpsiResult(const TH1& hminv) :
 282 AliAnalysisMuMuResult("",""),
 283 fNofRuns(1),
 284 fNofTriggers(-1),
 285 fMinv(0x0),
 286 fBin(),
 287 fRebin(0),
 288 fTriggerClass(),
 289 fEventSelection(),
 290 fPairSelection(),
 291 fCentralitySelection()
 292 {
 293   SetMinv(hminv);
 294 }
 295
 296 //_____________________________________________________________________________
 297 AliAnalysisMuMuJpsiResult::AliAnalysisMuMuJpsiResult(const TH1& hminv,
 298                                              const char* fitType,
 299                                              Int_t nrebin)
 300 :
 301 AliAnalysisMuMuResult(Form("%s:%d",fitType,nrebin),""),
 302 fNofRuns(1),
 303 fNofTriggers(-1),
 304 fMinv(0x0),
 305 fBin(),
 306 fRebin(nrebin),
 307 fTriggerClass(),
 308 fEventSelection(),
 309 fPairSelection(),
 310 fCentralitySelection()
 311 {
 312   SetMinv(hminv);
 313 }
 314
 315 //_____________________________________________________________________________
 316 AliAnalysisMuMuJpsiResult::AliAnalysisMuMuJpsiResult(const TH1& hminv,
 317                                              const char* triggerName,
 318                                              const char* eventSelection,
 319                                              const char* pairSelection,
 320                                              const char* centSelection,
 321                                              const AliAnalysisMuMuBinning::Range& bin)
 322 :
 323 AliAnalysisMuMuResult(Form("%s-%s-%s-%s",triggerName,eventSelection,pairSelection,centSelection),""),
 324 fNofRuns(1),
 325 fNofTriggers(-1),
 326 fMinv(0x0),
 327 fBin(bin),
 328 fRebin(1),
 329 fTriggerClass(triggerName),
 330 fEventSelection(eventSelection),
 331 fPairSelection(pairSelection),
 332 fCentralitySelection(centSelection)
 333 {
 334   SetMinv(hminv);
 335 }
 336
 337 //_____________________________________________________________________________
 338 AliAnalysisMuMuJpsiResult::AliAnalysisMuMuJpsiResult(const AliAnalysisMuMuJpsiResult& rhs)
 339 :
 340 AliAnalysisMuMuResult(rhs),
 341 fNofRuns(rhs.NofRuns()),
 342 fNofTriggers(rhs.NofTriggers()),
 343 fMinv(0x0),
 344 fBin(rhs.Bin()),
 345 fRebin(rhs.fRebin),
 346 fTriggerClass(rhs.fTriggerClass),
 347 fEventSelection(rhs.fEventSelection),
 348 fPairSelection(rhs.fPairSelection),
 349 fCentralitySelection(rhs.fCentralitySelection)
 350 {
 351   /// copy ctor
 352   /// Note that the mother is lost
 353   /// fKeys remains 0x0 so it will be recomputed if need be
 354
 355   if ( rhs.fMinv )
 356   {
 357     fMinv = static_cast<TH1*>(rhs.fMinv->Clone());
 358   }
 359 }
 360
 361 //_____________________________________________________________________________
 362 AliAnalysisMuMuJpsiResult& AliAnalysisMuMuJpsiResult::operator=(const AliAnalysisMuMuJpsiResult& rhs)
 363 {
 364   /// Assignment operator
 365
 366   if (this!=&rhs)
 367   {
 368     static_cast<AliAnalysisMuMuResult&>(*this) = static_cast<const AliAnalysisMuMuResult&>(rhs);
 369     delete fMinv;
 370
 371     if ( rhs.fMinv )
 372     {
 373       fMinv = static_cast<TH1*>(rhs.fMinv->Clone());
 374     }
 375
 376     fNofRuns = rhs.NofRuns();
 377     fNofTriggers = rhs.NofTriggers();
 378     fBin = rhs.Bin();
 379     fRebin = rhs.fRebin;
 380   }
 381
 382   return *this;
 383 }
 384
 385 //_____________________________________________________________________________
 386 AliAnalysisMuMuJpsiResult::~AliAnalysisMuMuJpsiResult()
 387 {
 388   // dtor
 389   delete fMinv;
 390 }
 391
 392 //_____________________________________________________________________________
 393 const AliAnalysisMuMuBinning::Range& AliAnalysisMuMuJpsiResult::Bin() const
 394 {
 395   /// Get the bin of this result
 396   if ( !Mother() ) return fBin;
 397   else return Mother()->Bin();
 398 }
 399
 400 //_____________________________________________________________________________
 401 TObject* AliAnalysisMuMuJpsiResult::Clone(const char* /*newname*/) const
 402 {
 403   /// Clone this result
 404   return new AliAnalysisMuMuJpsiResult(*this);
 405 }
 406
 407 //_____________________________________________________________________________
 408 Bool_t AliAnalysisMuMuJpsiResult::Correct(const AliAnalysisMuMuJpsiResult& other, const char* particle, const char* subResultName)
 409 {
 410   /// Assuming other has an AccxEff entry, correct this value by AccxEff of other
 411
 412   if ( HasValue(Form("Nof%s",particle)) )
 413   {
 414     if (!other.HasValue(Form("AccEff%s",particle),subResultName))
 415     {
 416       AliError(Form("Cannot correct as I do not find the AccEff%s value (subResultName=%s)!",particle,subResultName));
 417       return kFALSE;
 418     }
 419
 420     Double_t acc = other.GetValue(Form("AccEff%s",particle),subResultName);
 421     Double_t value = 0.0;
 422
 423     if ( acc > 0 ) value = GetValue(Form("Nof%s",particle)) / acc;
 424
 425     Double_t error = ErrorAB( GetValue(Form("Nof%s",particle)),
 426                               GetErrorStat(Form("Nof%s",particle)),
 427                               other.GetValue(Form("AccEff%s",particle),subResultName),
 428                               other.GetErrorStat(Form("AccEff%s",particle),subResultName) );
 429
 430     Set(Form("CorrNof%s",particle),value,error*value);
 431
 432     return kTRUE;
 433   }
 434   else
 435   {
 436     AliError(Form("Result does not have Nof%s : cannot correct it !",particle));
 437   }
 438   return kFALSE;
 439 }
 440
 441 //_____________________________________________________________________________
 442 Double_t AliAnalysisMuMuJpsiResult::CountParticle(const TH1& hminv, const char* particle, Double_t sigma)
 443 {
 444   /// Count the number of entries in an invariant mass range
 445
 446   const std::map<std::string,Double_t>& m = MassMap();
 447
 448   std::map<std::string,Double_t>::const_iterator it = m.find(particle);
 449
 450   if ( it == m.end() )
 451   {
 452     AliErrorClass(Form("Don't know the mass of particle %s",particle));
 453     return 0.0;
 454   }
 455
 456   Double_t mass = it->second;
 457
 458   if ( sigma < 0 )
 459   {
 460     AliDebugClass(1,Form("Oups. Got a sigma of %e for particle %s !",sigma,particle));
 461     return hminv.Integral();
 462   }
 463
 464   TAxis* x = hminv.GetXaxis();
 465
 466   Int_t b1 = x->FindBin(mass-sigma);
 467   Int_t b2 = x->FindBin(mass+sigma);
 468
 469   AliDebugClass(1,Form("hminv getentries %e integral %e",hminv.GetEntries(),hminv.Integral(b1,b2)));
 470
 471   return hminv.Integral(b1,b2);
 472 }
 473
 474 //_____________________________________________________________________________
 475 AliAnalysisMuMuJpsiResult* AliAnalysisMuMuJpsiResult::FitJpsiGCBE(TH1& h)
 476 {
 477   /// Fit Jpsi spectra with crystal ball + gaussian + exponential
 478
 479   std::cout << "Fit with jpsi alone (gaus + CB + expo)" << std::endl;
 480
 481   Int_t nrebin = fMinv->GetXaxis()->GetNbins() / h.GetXaxis()->GetNbins();
 482
 483   AliAnalysisMuMuJpsiResult* r = new AliAnalysisMuMuJpsiResult(h,"JPSIGCBE",nrebin);
 484
 485   TH1* hfit = r->Minv();
 486
 487   const Double_t xmin(1.0);
 488   const Double_t xmax(8.0);
 489
 490   TF1* fitTotal = new TF1("fitTotal",funcJpsiGCBE,xmin,xmax,10);
 491   fitTotal->SetParNames("cste","x0","sigma0","N","alpha","n","mean","sigma","expocste","exposlope");
 492
 493   fitTotal->SetParLimits(3,0,h.GetMaximum()*2); // N
 494
 495   const Double_t cbalpha(0.98);
 496   const Double_t cbn(5.2);
 497
 498   fitTotal->FixParameter(4,cbalpha);
 499   fitTotal->FixParameter(5,cbn);
 500
 501   fitTotal->SetParLimits(6,2.8,3.2); // mean
 502   fitTotal->SetParLimits(7,0.02,0.3); // sigma
 503
 504   TF1* fg = new TF1("fg","gaus",xmin,xmax);
 505
 506   hfit->Fit(fg,"","",0.75,3.0);
 507
 508   fitTotal->SetParameter(0,fg->GetParameter(0));
 509   fitTotal->SetParameter(1,fg->GetParameter(1));
 510   fitTotal->SetParameter(2,fg->GetParameter(2));
 511
 512   TF1* fexpo = new TF1("expo","expo",xmin,xmax);
 513
 514   hfit->Fit(fexpo,"","",3.5,5);
 515
 516   fitTotal->SetParameter(8,fexpo->GetParameter(0));
 517   fitTotal->SetParameter(9,fexpo->GetParameter(1));
 518
 519   fitTotal->SetParameter(3,h.GetMaximum()),
 520   fitTotal->SetParameter(4,cbalpha);
 521   fitTotal->SetParameter(5,cbn);
 522   fitTotal->SetParameter(6,3.15);
 523   fitTotal->SetParameter(7,0.1);
 524
 525   const char* fitOption = "QSI+";
 526
 527   TFitResultPtr fitResult = hfit->Fit(fitTotal,fitOption,"",2,5);
 528
 529   r->Set("MeanJpsi",fitTotal->GetParameter(6),fitTotal->GetParError(6));
 530   r->Set("SigmaJpsi",fitTotal->GetParameter(7),fitTotal->GetParError(7));
 531
 532   double m = r->GetValue("MeanJpsi");
 533   double s = r->GetValue("SigmaJpsi");
 534   double n = 3.0;
 535
 536   TF1* fcb = new TF1("fcb",funcCB,xmin,xmax,5);
 537   fcb->SetParameters(fitTotal->GetParameter(3),
 538                      fitTotal->GetParameter(4),
 539                      fitTotal->GetParameter(5),
 540                      fitTotal->GetParameter(6),
 541                      fitTotal->GetParameter(7));
 542
 543   fcb->SetLineColor(6);
 544   fcb->SetNpx(100);
 545   TLine* l1 = new TLine(m-n*s,0,m-n*s,fitTotal->GetParameter(3));
 546   TLine* l2 = new TLine(m+n*s,0,m+n*s,fitTotal->GetParameter(3));
 547   l1->SetLineColor(6);
 548   l2->SetLineColor(6);
 549   h.GetListOfFunctions()->Add(fcb);
 550   h.GetListOfFunctions()->Add(l1);
 551   h.GetListOfFunctions()->Add(l2);
 552
 553
 554   Double_t cbParameters[5];
 555   Double_t covarianceMatrix[5][5];
 556
 557   cbParameters[0] = fitTotal->GetParameter(3);
 558   cbParameters[1] = fitTotal->GetParameter(4);
 559   cbParameters[2] = fitTotal->GetParameter(5);
 560   cbParameters[3] = fitTotal->GetParameter(6);
 561   cbParameters[4] = fitTotal->GetParameter(7);
 562
 563   for ( int iy = 0; iy < 5; ++iy )
 564   {
 565     for ( int ix = 0; ix < 5; ++ix )
 566     {
 567       covarianceMatrix[ix][iy] = (fitResult->GetCovarianceMatrix())(ix+3,iy+3);
 568     }
 569   }
 570
 571   double njpsi = fcb->Integral(m-n*s,m+n*s)/h.GetBinWidth(1);
 572
 573   double nerr = fcb->IntegralError(m-n*s,m+n*s,&cbParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(1);
 574
 575   r->Set("NofJpsi",njpsi,nerr);
 576
 577   return r;
 578 }
 579
 580 //_____________________________________________________________________________
 581 AliAnalysisMuMuJpsiResult* AliAnalysisMuMuJpsiResult::FitJpsi(TH1& h)
 582 {
 583   /// Fit Jpsi spectra using extended crystall ball (CB2) with free tails
 584
 585   StdoutToAliDebug(1,std::cout << "Fit with jpsi alone" << std::endl;);
 586
 587   Int_t nrebin = fMinv->GetXaxis()->GetNbins() / h.GetXaxis()->GetNbins();
 588
 589   AliAnalysisMuMuJpsiResult* r = new AliAnalysisMuMuJpsiResult(h,"JPSI",nrebin);
 590
 591   TH1* hfit = r->Minv();
 592
 593   const Double_t xmin(1.5);
 594   const Double_t xmax(8.0);
 595
 596   TF1* fitTotal = new TF1("fitTotal",funcCB2,xmin,xmax,7);
 597   fitTotal->SetParNames("N","alphaLow","nLow","mean","sigma","alphaUp","nUp");
 598   fitTotal->SetParameters(h.GetMaximum(),1,5,3.1,0.07,1.5,3);
 599   fitTotal->SetParLimits(0,0,h.GetMaximum()*2); // N
 600   fitTotal->SetParLimits(1,0,10); // alpha
 601   fitTotal->SetParLimits(2,0.1,10); // n
 602   fitTotal->SetParLimits(3,3,3.15); // mean
 603   fitTotal->SetParLimits(4,0.01,1); // sigma
 604   fitTotal->SetParLimits(5,0,10); // alpha
 605   fitTotal->SetParLimits(6,0.1,10); // n
 606
 607   hfit->Fit(fitTotal,"QSER+","",2,5);
 608
 609
 610   r->Set("MeanJpsi",fitTotal->GetParameter(3),fitTotal->GetParError(3));
 611   r->Set("SigmaJpsi",fitTotal->GetParameter(4),fitTotal->GetParError(4));
 612
 613   double m = r->GetValue("MeanJpsi");
 614   double s = r->GetValue("SigmaJpsi");
 615   double n = 10.0;
 616
 617   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));
 618
 619   return r;
 620 }
 621
 622 //_____________________________________________________________________________
 623 AliAnalysisMuMuJpsiResult* AliAnalysisMuMuJpsiResult::FitJpsiCB2VWG(const TH1& h)
 624 {
 625   /// Fit Jpsi spectra using extended crystal ball (CB2) + variable width gaussian (VWG)
 626
 627   StdoutToAliDebug(1,std::cout << "Fit with jpsi VWG" << std::endl;);
 628
 629   Int_t nrebin = fMinv->GetXaxis()->GetNbins() / h.GetXaxis()->GetNbins();
 630
 631   AliAnalysisMuMuJpsiResult* r = new AliAnalysisMuMuJpsiResult(h,"JPSICB2VWG",nrebin);
 632
 633
 634   TH1* hfit = r->Minv();
 635
 636   const Double_t xmin(2.0);
 637   const Double_t xmax(5.0);
 638
 639 //  // gaussian variable width
 640 //  Double_t sigma = par[2]+par[3]*((x[0]-par[1])/par[1]);
 641 //  return par[0]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/(2.*sigma*sigma));
 642 //  Double_t CrystalBallExtended(Double_t *x,Double_t *par)
 643 //  //par[0] = Normalization
 644 //  //par[1] = mean
 645 //  //par[2] = sigma
 646 //  //par[3] = alpha
 647 //  //par[4] = n
 648 //  //par[5] = alpha'
 649 //  //par[6] = n'
 650
 651   TF1* fitTotal = new TF1("fitTotal",fitFunctionCB2VWG,xmin,xmax,11);
 652   fitTotal->SetParNames("kVWG","mVWG","sVWG1","sVWG2","norm","mean","sigma","alpha","n","alpha'","n'");
 653
 654   fitTotal->SetParameter(0, 10000.); // kVWG
 655   fitTotal->SetParameter(1, 1.9); // mVWG
 656
 657   fitTotal->SetParameter(2, 0.5); // sVWG1
 658   fitTotal->SetParLimits(2, 0., 100.);
 659
 660   fitTotal->SetParameter(3, 0.3); // sVWG2
 661   fitTotal->SetParLimits(3, 0., 100.);
 662
 663   fitTotal->SetParameter(4, h.GetMaximum()); // norm
 664
 665   fitTotal->SetParameter(5, 3.1); // mean
 666   fitTotal->SetParLimits(5, 3.0, 3.2);
 667
 668   fitTotal->SetParameter(6, 0.08); // sigma
 669   fitTotal->SetParLimits(6, 0.04, 0.20);
 670
 671   fitTotal->SetParameter(7,1.0); // alpha
 672   fitTotal->SetParameter(8,5); // n
 673   fitTotal->SetParameter(9,2.0); // alpha'
 674   fitTotal->SetParameter(10,4); // n'
 675
 676 //  fitTotal->FixParameter(7, 0.93);
 677 //  fitTotal->FixParameter(8, 5.59);
 678 //  fitTotal->FixParameter(9, 2.32);
 679 //  fitTotal->FixParameter(10, 3.39);
 680 //  fitTotal->SetParameter(11, 10.);
 681
 682   const char* fitOption = "QSIER"; //+";
 683
 684   TFitResultPtr fitResult = hfit->Fit(fitTotal,fitOption,"");
 685
 686   r->Set("MeanJpsi",fitTotal->GetParameter(5),fitTotal->GetParError(5));
 687   r->Set("SigmaJpsi",fitTotal->GetParameter(6),fitTotal->GetParError(6));
 688
 689   double m = r->GetValue("MeanJpsi");
 690   double s = r->GetValue("SigmaJpsi");
 691   double n = 3.0;
 692
 693   TF1* fcb = new TF1("fcb",CrystalBallExtended,xmin,xmax,7);
 694   fcb->SetParameters(fitTotal->GetParameter(4),
 695                      fitTotal->GetParameter(5),
 696                      fitTotal->GetParameter(6),
 697                      fitTotal->GetParameter(7),
 698                      fitTotal->GetParameter(8),
 699                      fitTotal->GetParameter(9),
 700                      fitTotal->GetParameter(10));
 701
 702
 703   fcb->SetLineColor(1);
 704   fcb->SetNpx(1000);
 705   TLine* l1 = new TLine(m-n*s,0,m-n*s,fitTotal->GetParameter(4));
 706   TLine* l2 = new TLine(m+n*s,0,m+n*s,fitTotal->GetParameter(4));
 707   l1->SetLineColor(6);
 708   l2->SetLineColor(6);
 709   hfit->GetListOfFunctions()->Add(fcb);
 710   hfit->GetListOfFunctions()->Add(l1);
 711   hfit->GetListOfFunctions()->Add(l2);
 712
 713   Double_t cbParameters[7];
 714   Double_t covarianceMatrix[7][7];
 715
 716   for ( int ix = 0; ix < 7; ++ix )
 717   {
 718     cbParameters[ix] = fitTotal->GetParameter(ix+4);
 719   }
 720
 721   for ( int iy = 0; iy < 5; ++iy )
 722   {
 723     for ( int ix = 0; ix < 5; ++ix )
 724     {
 725       covarianceMatrix[ix][iy] = (fitResult->GetCovarianceMatrix())(ix+4,iy+4);
 726     }
 727   }
 728
 729   double njpsi = fcb->Integral(m-n*s,m+n*s)/h.GetBinWidth(1);
 730   double nerr = fcb->IntegralError(m-n*s,m+n*s,&cbParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(1);
 731
 732   r->Set("NofJpsi",njpsi,nerr);
 733
 734   return r;
 735 }
 736
 737 //_____________________________________________________________________________
 738 AliAnalysisMuMuJpsiResult* AliAnalysisMuMuJpsiResult::FitJpsi2CB2VWG(const TH1& h,
 739                                                              Double_t alphaLow,
 740                                                              Double_t nLow,
 741                                                              Double_t alphaUp,
 742                                                              Double_t nUp)
 743 {
 744   /// Fit using extended crystal ball + variable width gaussian
 745
 746   StdoutToAliDebug(1,std::cout << Form("Fit with jpsi + psiprime VWG alphaLow=%5.2f nLow=%5.2f alphaUp=%5.2f nUp=%5.2f",
 747                                        alphaLow,nLow,alphaUp,nUp) << std::endl;);
 748
 749   Int_t nrebin = fMinv->GetXaxis()->GetNbins() / h.GetXaxis()->GetNbins();
 750
 751   TString resultName("JPSI2CB2VWG");
 752   if ( alphaLow > 0 )
 753   {
 754     resultName += TString::Format("alphaLow=%5.2f",alphaLow);
 755   }
 756   if ( nLow > 0 )
 757   {
 758     resultName += TString::Format("nLow=%5.2f",nLow);
 759   }
 760   if ( alphaUp > 0 )
 761   {
 762     resultName += TString::Format("alphaUp=%5.2f",alphaUp);
 763   }
 764   if ( nUp > 0 )
 765   {
 766     resultName += TString::Format("nUp=%5.2f",nUp);
 767   }
 768   resultName.ReplaceAll(" ","");
 769
 770   AliAnalysisMuMuJpsiResult* r = new AliAnalysisMuMuJpsiResult(h,resultName.Data(),nrebin);
 771
 772   TH1* hfit = r->Minv();
 773
 774   const Double_t xmin(2.2);
 775   const Double_t xmax(5.0);
 776
 777   TF1* fitTotal = new TF1("fitTotal",func2CB2VWG,xmin,xmax,12);
 778   fitTotal->SetParNames("kVWG","mVWG","sVWG1","sVWG2","kPsi","mPsi","sPsi","alPsi","nlPsi","auPsi","nuPsi");
 779   fitTotal->SetParName(11, "kPsi'");
 780
 781   fitTotal->SetParameter(0, 10000.);
 782   fitTotal->SetParameter(1, 1.9);
 783   fitTotal->SetParameter(2, 0.5);
 784   fitTotal->SetParLimits(2, 0., 100.);
 785   fitTotal->SetParameter(3, 0.3);
 786   fitTotal->SetParLimits(3, 0., 100.);
 787   fitTotal->SetParameter(4, 100.);
 788   fitTotal->SetParameter(5, 3.1);
 789   fitTotal->SetParLimits(5, 3.08, 3.2);
 790   fitTotal->SetParameter(6, 0.08);
 791   fitTotal->SetParLimits(6, 0.05, 0.15);
 792
 793 //  r = FitJpsi2CB2VWG(*hminv,0.93,5.59,2.32,3.39);
 794
 795   if ( alphaLow > 0 )
 796   {
 797     fitTotal->FixParameter(7, alphaLow);
 798   }
 799   else
 800   {
 801     fitTotal->SetParameter(7,0.9);
 802     fitTotal->SetParLimits(7,0.1,10.0);
 803   }
 804
 805   if ( nLow > 0 )
 806   {
 807     fitTotal->FixParameter(8, nLow);
 808   }
 809   else
 810   {
 811     fitTotal->SetParameter(8,5.0);
 812     fitTotal->SetParLimits(8,0.0,10.0);
 813   }
 814
 815   if ( alphaUp > 0 )
 816   {
 817     fitTotal->FixParameter(9, alphaUp);
 818   }
 819   else
 820   {
 821     fitTotal->SetParameter(9, 2.0);
 822     fitTotal->SetParLimits(9,0.1,10.0);
 823   }
 824
 825   if ( nUp > 0 )
 826   {
 827     fitTotal->FixParameter(10, nUp);
 828   }
 829   else
 830   {
 831     fitTotal->SetParameter(10,3.0);
 832     fitTotal->SetParLimits(10,0.0,10.0);
 833   }
 834
 835   fitTotal->SetParameter(11, 10.);
 836
 837   const char* fitOption = "QSER"; //+";
 838
 839   TFitResultPtr fitResult = hfit->Fit(fitTotal,fitOption,"");
 840
 841   r->Set("MeanJpsi",fitTotal->GetParameter(5),fitTotal->GetParError(5));
 842   r->Set("SigmaJpsi",fitTotal->GetParameter(6),fitTotal->GetParError(6));
 843
 844   double m = r->GetValue("MeanJpsi");
 845   double s = r->GetValue("SigmaJpsi");
 846   double n = 3.0;
 847
 848   TF1* fcb = new TF1("fcb",CrystalBallExtended,xmin,xmax,7);
 849   fcb->SetParameters(fitTotal->GetParameter(4),
 850                      fitTotal->GetParameter(5),
 851                      fitTotal->GetParameter(6),
 852                      fitTotal->GetParameter(7),
 853                      fitTotal->GetParameter(8),
 854                      fitTotal->GetParameter(9),
 855                      fitTotal->GetParameter(10));
 856
 857
 858   fcb->SetLineColor(1);
 859   fcb->SetNpx(1000);
 860   TLine* l1 = new TLine(m-n*s,0,m-n*s,fitTotal->GetParameter(4));
 861   TLine* l2 = new TLine(m+n*s,0,m+n*s,fitTotal->GetParameter(4));
 862   l1->SetLineColor(6);
 863   l2->SetLineColor(6);
 864   hfit->GetListOfFunctions()->Add(fcb);
 865   hfit->GetListOfFunctions()->Add(l1);
 866   hfit->GetListOfFunctions()->Add(l2);
 867
 868   Double_t cbParameters[7];
 869   Double_t covarianceMatrix[7][7];
 870
 871   for ( int ix = 0; ix < 7; ++ix )
 872   {
 873     cbParameters[ix] = fitTotal->GetParameter(ix+4);
 874   }
 875
 876   for ( int iy = 0; iy < 5; ++iy )
 877   {
 878     for ( int ix = 0; ix < 5; ++ix )
 879     {
 880       covarianceMatrix[ix][iy] = (fitResult->GetCovarianceMatrix())(ix+4,iy+4);
 881     }
 882   }
 883
 884   double njpsi = fcb->Integral(m-n*s,m+n*s)/h.GetBinWidth(1);
 885   double nerr = fcb->IntegralError(m-n*s,m+n*s,&cbParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(1);
 886
 887   r->Set("NofJpsi",njpsi,nerr);
 888
 889   return r;
 890 }
 891
 892 //_____________________________________________________________________________
 893 AliAnalysisMuMuJpsiResult* AliAnalysisMuMuJpsiResult::FitJpsiNA48(const TH1& h)
 894 {
 895   /// fit using functional form from Ruben Shahoyan's thesis (2001) (eq. 4.8.)
 896
 897   StdoutToAliDebug(1,std::cout << "Fit with jpsi NA50 Ruben eq. 4.8" << std::endl;);
 898
 899   Int_t nrebin = fMinv->GetXaxis()->GetNbins() / h.GetXaxis()->GetNbins();
 900
 901   AliAnalysisMuMuJpsiResult* r = new AliAnalysisMuMuJpsiResult(h,"JPSINA",nrebin);
 902
 903   TH1* hfit = r->Minv();
 904
 905   const Double_t xmin(2.0);
 906   const Double_t xmax(5.0);
 907
 908   TF1* fitTotal = new TF1("fitTotal",funcJpsiNA48,xmin,xmax,12);
 909
 910   fitTotal->SetParName( 0, "c1");
 911   fitTotal->FixParameter(0,0.97);
 912
 913   fitTotal->SetParName( 1, "c2");
 914   fitTotal->FixParameter(1,1.05);
 915
 916   fitTotal->SetParName( 2, "d1");
 917   fitTotal->SetParameter(2,0.0);
 918   fitTotal->SetParLimits(2,0,1);
 919
 920   fitTotal->SetParName( 3, "d2");
 921   fitTotal->SetParameter(3,0.0);
 922   fitTotal->SetParLimits(3,0,1);
 923
 924   fitTotal->SetParName( 4, "g1");
 925   fitTotal->SetParameter(4,0.0);
 926   fitTotal->SetParLimits(4,0.0,1);
 927
 928   fitTotal->SetParName( 5, "g2");
 929   fitTotal->SetParameter(5,0.0);
 930   fitTotal->SetParLimits(5,0.0,1);
 931
 932   fitTotal->SetParName( 6, "m0");
 933   fitTotal->SetParameter(6,3.1);
 934   fitTotal->SetParLimits(6,2.8,3.4);
 935
 936   fitTotal->SetParName( 7, "sigma1");
 937   fitTotal->SetParameter(7,0.05);
 938   fitTotal->SetParLimits(7,0.01,0.2);
 939
 940   fitTotal->SetParName( 8, "sigma2");
 941   fitTotal->SetParameter(8,0.05);
 942   fitTotal->SetParLimits(8,0.01,0.2);
 943
 944   fitTotal->SetParName( 9, "b1");
 945   fitTotal->SetParameter(9,1.0);
 946   fitTotal->SetParLimits(9,0,1);
 947
 948   fitTotal->SetParName(10, "b2");
 949   fitTotal->SetParameter(10,1.0);
 950   fitTotal->SetParLimits(10,0,1);
 951
 952   fitTotal->SetParName(11, "norm");
 953   fitTotal->SetParameter(11,h.GetMaximum());
 954
 955   const char* fitOption = "QSIER"; //+";
 956
 957   TFitResultPtr fitResult = hfit->Fit(fitTotal,fitOption,"");
 958
 959   r->Set("MeanJpsi",fitTotal->GetParameter(6),fitTotal->GetParError(6));
 960   r->Set("SigmaJpsi",
 961          fitTotal->GetParameter(7)+fitTotal->GetParameter(8),
 962          0.0);
 963
 964   double m = r->GetValue("MeanJpsi");
 965   double s = r->GetValue("SigmaJpsi");
 966   double n = 3.0;
 967
 968   TLine* l1 = new TLine(m-n*s,0,m-n*s,fitTotal->GetParameter(11));
 969   TLine* l2 = new TLine(m+n*s,0,m+n*s,fitTotal->GetParameter(11));
 970   l1->SetLineColor(6);
 971   l2->SetLineColor(6);
 972   hfit->GetListOfFunctions()->Add(l1);
 973   hfit->GetListOfFunctions()->Add(l2);
 974
 975   double njpsi = fitTotal->Integral(m-n*s,m+n*s)/h.GetBinWidth(1);
 976   double nerr = fitTotal->IntegralError(m-n*s,m+n*s)/h.GetBinWidth(1);
 977
 978   r->Set("NofJpsi",njpsi,nerr);
 979
 980   return r;
 981 }
 982
 983 //_____________________________________________________________________________
 984 Bool_t AliAnalysisMuMuJpsiResult::AddFit(const char* fitType, Int_t npar, Double_t* par)
 985 {
 986   // Add a fit to this result
 987
 988   TString msg(Form("fitType=%s npar=%d par[]=",fitType,npar));
 989
 990   for ( Int_t i = 0; i < npar; ++i )
 991   {
 992     msg += TString::Format("%e,",par[i]);
 993   }
 994
 995   msg += TString::Format(" minv=%p %d",fMinv,fMinv?TMath::Nint(fMinv->GetEntries()):0);
 996
 997   if ( !fMinv ) return kFALSE;
 998
 999   TString sFitType(fitType);
1000   sFitType.ToUpper();
1001   Int_t nrebin(1);
1002
1003   if (sFitType.CountChar(':'))
1004   {
1005     Int_t index = sFitType.Index(":");
1006     nrebin = TString(sFitType(index+1,sFitType.Length()-index-1)).Atoi();
1007     sFitType = sFitType(0,index);
1008   }
1009
1010   msg += TString::Format(" nrebin=%d",nrebin);
1011
1012   AliDebug(1,msg.Data());
1013
1014
1015   if ( fMinv->GetEntries()<100 && !sFitType.Contains("COUNT")) return kFALSE;
1016
1017   TH1* hminv = static_cast<TH1*>(fMinv->Clone());
1018
1019   hminv->Rebin(nrebin);
1020   hminv->SetDirectory(0);
1021
1022   AliAnalysisMuMuJpsiResult* r(0x0);
1023
1024   if ( sFitType=="PSI1")
1025   {
1026     r = FitJpsi(*hminv);
1027   }
1028   else if ( sFitType == "PSILOW")
1029   {
1030     r = FitJpsi2CB2VWG(*hminv,-1,-1,-1,-1); // free tails
1031   }
1032   else if ( sFitType == "PSILOWMCTAILS" )
1033   {
1034     if ( npar!= 4 )
1035     {
1036       AliError("Cannot use PSILOWMCTAILS without being given the MC tails !");
1037       delete hminv;
1038       return kFALSE;
1039     }
1040     r = FitJpsi2CB2VWG(*hminv,par[0],par[1],par[2],par[3]);
1041     if (r)
1042     {
1043       r->SetAlias("MCTAILS");
1044     }
1045   }
1046   else if ( sFitType.BeginsWith("PSILOWALPHA") )
1047   {
1048     Float_t lpar[] = { 0.0, 0.0, 0.0, 0.0 };
1049
1050     AliDebug(1,Form("sFitType=%s",sFitType.Data()));
1051
1052     sscanf(sFitType.Data(),"PSILOWALPHALOW%fNLOW%fALPHAUP%fNUP%f",
1053            &lpar[0],&lpar[1],&lpar[2],&lpar[3]);
1054
1055     AliDebug(1,Form("PSILOW ALPHALOW=%f NLOW=%f ALPHAUP=%f NUP=%f",lpar[0],lpar[1],lpar[2],lpar[3]));
1056
1057     if ( lpar[0] == 0.0 && lpar[1] == 0.0 && lpar[0] == 0.0 && lpar[1] == 0.0 )
1058     {
1059       AliError("Cannot work with zero tails !");
1060     }
1061     else
1062     {
1063       r = FitJpsi2CB2VWG(*hminv,lpar[0],lpar[1],lpar[2],lpar[3]);
1064     }
1065   }
1066   else if ( sFitType == "COUNTJPSI" )
1067   {
1068     r = new AliAnalysisMuMuJpsiResult(*hminv,"COUNTJPSI",1);
1069     Double_t n = CountParticle(*hminv,"Jpsi");
1070     r->Set("NofJpsi",n,TMath::Sqrt(n));
1071   }
1072
1073
1074   if ( r )
1075   {
1076     StdoutToAliDebug(1,r->Print(););
1077     r->SetBin(Bin());
1078     r->SetNofTriggers(NofTriggers());
1079     r->SetNofRuns(NofRuns());
1080
1081     AdoptSubResult(r);
1082   }
1083
1084   delete hminv;
1085
1086   return (r!=0x0);
1087 }
1088
1089 //_____________________________________________________________________________
1090 Long64_t AliAnalysisMuMuJpsiResult::Merge(TCollection* list)
1091 {
1092   /// Merge method
1093   ///
1094   /// Merge a list of AliAnalysisMuMuJpsiResult objects with this
1095   /// Returns the number of merged objects (including this).
1096   ///
1097   /// Note that the merging is to be understood here as an weighted mean operation
1098   ///
1099   /// FIXME ! (compared to base class Merge, should only Minv merging ?)
1100
1101   AliError("Implement me !");
1102   if (!list) return 0;
1103
1104   if (list->IsEmpty()) return 1;
1105
1106   return 0;
1107
1108 }
1109
1110 //_____________________________________________________________________________
1111 Int_t AliAnalysisMuMuJpsiResult::NofRuns() const
1112 {
1113   /// Get the number of runs
1114   if ( !Mother() ) return fNofRuns;
1115   else return Mother()->NofRuns();
1116 }
1117
1118 //_____________________________________________________________________________
1119 Int_t AliAnalysisMuMuJpsiResult::NofTriggers() const
1120 {
1121   /// Get the number of triggers
1122
1123   if ( !Mother() ) return fNofTriggers;
1124   else return Mother()->NofTriggers();
1125 }
1126
1127 //_____________________________________________________________________________
1128 void AliAnalysisMuMuJpsiResult::Print(Option_t* opt) const
1129 {
1130   /// printout
1131
1132   std::cout << Form("NRUNS %d - NTRIGGER %10d - %s",
1133                                       NofRuns(),
1134                                       NofTriggers(),
1135                                       fBin.AsString().Data());
1136
1137   AliAnalysisMuMuResult::Print(opt);
1138 }
1139
1140 //_____________________________________________________________________________
1141 void AliAnalysisMuMuJpsiResult::PrintValue(const char* key, const char* opt, Double_t value, Double_t errorStat,
1142                                            Double_t rms) const
1143 {
1144   /// exclude the particles with zero stat
1145
1146   const std::map<std::string,Double_t>& m = MassMap();
1147
1148   for( std::map<std::string,Double_t>::const_iterator it = m.begin(); it != m.end(); ++it )
1149   {
1150     TString particle(it->first.c_str());
1151
1152     if (TString(key).Contains(particle.Data()))
1153     {
1154       if ( GetValue("Nof%s",particle.Data()) <= 0.0 ) return;
1155     }
1156   }
1157
1158   AliAnalysisMuMuResult::PrintValue(key,opt,value,errorStat,rms);
1159
1160 }
1161
1162 //_____________________________________________________________________________
1163 void AliAnalysisMuMuJpsiResult::PrintParticle(const char* particle, const char* opt) const
1164 {
1165   /// Print all information about one particule type
1166
1167   Double_t npart = GetValue(Form("Nof%s",particle));
1168   if (npart<=0) return;
1169
1170
1171   std::cout << opt << Form("\t%s",particle) << std::endl;
1172
1173   //  Double_t npartError = GetErrorStat(Form("Nof%s",particle));
1174 //  std::cout << opt << Form("\t\t%20s %9.2f +- %5.2f","Count",npart,npartError) << std::endl;
1175
1176   TIter next(Keys());
1177   TObjString* key;
1178
1179   while ( ( key = static_cast<TObjString*>(next()) ) )
1180   {
1181     if ( !key->String().Contains(particle) ) continue;
1182
1183     PrintValue(key->String(),opt,GetValue(key->String()),GetErrorStat(key->String()),GetRMS(key->String()));
1184   }
1185 }
1186
1187 //_____________________________________________________________________________
1188 void AliAnalysisMuMuJpsiResult::SetBin(const AliAnalysisMuMuBinning::Range& bin)
1189 {
1190   /// Set the bin
1191
1192   if (!Mother()) fBin = bin;
1193   else Mother()->SetBin(bin);
1194 }
1195
1196 //_____________________________________________________________________________
1197 void AliAnalysisMuMuJpsiResult::SetNofInputParticles(const TH1& hminv)
1198 {
1199   /// Set the number of input particle from the invariant mass spectra
1200
1201   static const char* particleNames[] = { "Jpsi" , "PsiPrime", "Upsilon","UpsilonPrime" };
1202
1203   const std::map<std::string,Double_t>& m = MassMap();
1204
1205   for ( Int_t i = 0; i < 4; ++i )
1206   {
1207     std::map<std::string,Double_t>::const_iterator it = m.find(particleNames[i]);
1208
1209     Double_t sigma(-1.0);
1210
1211     if (it != m.end() )
1212     {
1213       sigma = it->second*0.1;
1214     }
1215
1216     Double_t n = CountParticle(hminv,particleNames[i],sigma);
1217
1218     AliDebug(1,Form("i=%d particle %s n %e",i,particleNames[i],n));
1219
1220     if ( n > 0 )
1221     {
1222       SetNofInputParticles(particleNames[i],TMath::Nint(n));
1223     }
1224   }
1225 }
1226
1227 //_____________________________________________________________________________
1228 void AliAnalysisMuMuJpsiResult::SetNofInputParticles(const char* particle, int n)
1229 {
1230   /// Set the number of input particles (so it is a MC result)
1231   /// and (re)compute the AccxEff values
1232
1233   Set(Form("NofInput%s",particle),n,TMath::Sqrt(n));
1234
1235   if (n<=0)
1236   {
1237     Set(Form("AccEff%s",particle),0,0);
1238     return;
1239   }
1240
1241   Double_t npart = GetValue(Form("Nof%s",particle));
1242   Double_t npartErr  = GetErrorStat(Form("Nof%s",particle));
1243   Double_t ninput = GetValue(Form("NofInput%s",particle));
1244   Double_t ninputErr = GetErrorStat(Form("NofInput%s",particle));
1245
1246   Set(Form("AccEff%s",particle),
1247       npart/ninput,
1248       (npart/ninput)*ErrorAB(npart,npartErr,ninput,ninputErr));
1249
1250   TIter next(SubResults());
1251   AliAnalysisMuMuJpsiResult* r;
1252
1253   while ( ( r = static_cast<AliAnalysisMuMuJpsiResult*>(next())) )
1254   {
1255     r->Set(Form("NofInput%s",particle),n,TMath::Sqrt(n));
1256
1257     npart = r->GetValue(Form("Nof%s",particle));
1258     npartErr = r->GetErrorStat(Form("Nof%s",particle));
1259
1260     r->Set(Form("AccEff%s",particle),
1261            npart/ninput,
1262            (npart/ninput)*ErrorAB(npart,npartErr,ninput,ninputErr));
1263
1264   }
1265 }
1266
1267 //_____________________________________________________________________________
1268 void AliAnalysisMuMuJpsiResult::SetNofRuns(Int_t n)
1269 {
1270   if ( !Mother() ) fNofRuns=n;
1271   else Mother()->SetNofRuns(n);
1272 }
1273
1274 //_____________________________________________________________________________
1275 void AliAnalysisMuMuJpsiResult::SetNofTriggers(Int_t n)
1276 {
1277   if ( !Mother() ) fNofTriggers=n;
1278   else Mother()->SetNofTriggers(n);
1279 }
1280
1281 //_____________________________________________________________________________
1282 void AliAnalysisMuMuJpsiResult::SetMinv(const TH1& hminv)
1283 {
1284     /// Set the inv. mass spectrum to be fitted.
1285   static UInt_t n(0);
1286
1287   delete fMinv;
1288   fMinv = static_cast<TH1*>(hminv.Clone(Form("Minv%u",n++)));
1289   fMinv->SetDirectory(0);
1290 }