PWG/muondep/AliAnalysisMuMu.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 "AliAnalysisMuMu.h"
  20 #include "AliAnalysisMuMuBase.h" //Just to have avaliable the MCInputPrefix() method
  21
  22 #include "AliAnalysisMuMuBinning.h"
  23 #include "AliAnalysisMuMuConfig.h"
  24 #include "AliAnalysisMuMuFnorm.h"
  25 #include "AliAnalysisMuMuGraphUtil.h"
  26 #include "AliAnalysisMuMuJpsiResult.h"
  27 #include "AliAnalysisMuMuSpectra.h"
  28 #include "AliAnalysisTriggerScalers.h"
  29 #include "AliCounterCollection.h"
  30 #include "AliHistogramCollection.h"
  31 #include "AliLog.h"
  32 #include "AliMergeableCollection.h"
  33 #include "Riostream.h"
  34 #include "TArrayL64.h"
  35 #include "TASImage.h"
  36 #include "TAxis.h"
  37 #include "TCanvas.h"
  38 #include "TColor.h"
  39 #include "TF1.h"
  40 #include "TFile.h"
  41 #include "TGraphErrors.h"
  42 #include "TGrid.h"
  43 #include "TH1.h"
  44 #include "TH2.h"
  45 #include "TProfile.h"
  46 #include "THashList.h"
  47 #include "TKey.h"
  48 #include "TLegend.h"
  49 #include "TLegendEntry.h"
  50 #include "TLine.h"
  51 #include "TList.h"
  52 #include "TMap.h"
  53 #include "TMath.h"
  54 #include "TObjArray.h"
  55 #include "TObjString.h"
  56 #include "TParameter.h"
  57 #include "TPaveText.h"
  58 #include "TPRegexp.h"
  59 #include "TROOT.h"
  60 #include "TStopwatch.h"
  61 #include "TStyle.h"
  62 #include "TSystem.h"
  63 #include <cassert>
  64 #include <map>
  65 #include <set>
  66 #include <string>
  67 #include "TLatex.h"
  68
  69 ClassImp(AliAnalysisMuMu)
  70
  71 //_____________________________________________________________________________
  72 AliAnalysisMuMu::AliAnalysisMuMu(const char* filename, const char* associatedSimFileName, const char* associatedSimFileName2, const char* beamYear) : TObject(),
  73 fFilename(filename),
  74 fCounterCollection(0x0),
  75 fBinning(0x0),
  76 fMergeableCollection(0x0),
  77 fRunNumbers(),
  78 fCorrectionPerRun(0x0),
  79 fAssociatedSimulation(0x0),
  80 fAssociatedSimulation2(0x0),
  81 fParticleName(""),
  82 fConfig(0x0)
  83 {
  84   // ctor
  85
  86   GetCollections(fFilename,fMergeableCollection,fCounterCollection,fBinning,fRunNumbers);
  87
  88   if ( IsSimulation() )
  89   {
  90     TString sFilename(filename);
  91     if ( sFilename.Contains("JPSI",TString::kIgnoreCase) ) SetParticleName("JPsi"); //FIXME: DO a method for this
  92     else if ( sFilename.Contains("PSIP",TString::kIgnoreCase) ) SetParticleName("PsiP");
  93     else AliError("Unknown Particle Name in simulation");
  94   }
  95
  96   if ( fCounterCollection )
  97   {
  98     if ( strlen(associatedSimFileName) )
  99     {
 100       fAssociatedSimulation = new AliAnalysisMuMu(associatedSimFileName);
 101
 102       TString sAssociatedSimFileName(associatedSimFileName);
 103       if ( sAssociatedSimFileName.Contains("JPSI",TString::kIgnoreCase) ) fAssociatedSimulation->SetParticleName("JPsi");
 104       else if ( sAssociatedSimFileName.Contains("PSIP",TString::kIgnoreCase) ) fAssociatedSimulation->SetParticleName("PsiP");
 105       else AliError("Unknown Particle Name in associated simulation");
 106     }
 107
 108     if ( strlen(associatedSimFileName2) )
 109     {
 110       fAssociatedSimulation2 = new AliAnalysisMuMu(associatedSimFileName2);
 111
 112       TString sAssociatedSimFileName2(associatedSimFileName2);
 113       if ( sAssociatedSimFileName2.Contains("JPSI",TString::kIgnoreCase) ) fAssociatedSimulation2->SetParticleName("JPsi");
 114       else if ( sAssociatedSimFileName2.Contains("PSIP",TString::kIgnoreCase) ) fAssociatedSimulation2->SetParticleName("PsiP");
 115       else AliError("Unknown Particle Name in associated simulation 2");
 116     }
 117
 118     fConfig = new AliAnalysisMuMuConfig(beamYear);
 119   }
 120 }
 121
 122 //_____________________________________________________________________________
 123 AliAnalysisMuMu::~AliAnalysisMuMu()
 124 {
 125   // dtor
 126
 127   if ( fAssociatedSimulation )
 128   {
 129     fAssociatedSimulation->Update();
 130   }
 131   if ( fAssociatedSimulation2 )
 132   {
 133     fAssociatedSimulation2->Update();
 134   }
 135
 136   Update();
 137
 138   delete fCounterCollection;
 139   delete fBinning;
 140   delete fMergeableCollection;
 141   delete fCorrectionPerRun;
 142   delete fAssociatedSimulation;
 143   delete fAssociatedSimulation2;
 144   delete fConfig;
 145 }
 146
 147 //_____________________________________________________________________________
 148 void AliAnalysisMuMu::BasicCounts(Bool_t detailTriggers,
 149                                   ULong64_t* totalNmb,
 150                                   ULong64_t* totalNmsl,
 151                                   ULong64_t* totalNmul)
 152 {
 153   // Report of some basic numbers, like number of MB and MUON triggers,
 154   // both before and after physics selection, and comparison with
 155   // the total number of such triggers (taken from the OCDB scalers)
 156   // if requested.
 157   //
 158   // filename is assumed to be a root filecontaining a list containing
 159   //    an AliCounterCollection (or directly an AliCounterCollection)
 160   //
 161   // if detailTriggers is kTRUE, each kind of (MB,MUL,MSL) is counted separately
 162   //
 163
 164   if (!fMergeableCollection || !fCounterCollection) return;
 165
 166   TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
 167   TIter nextRun(runs);
 168
 169   TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
 170   TIter nextTrigger(triggers);
 171
 172   TObjArray* events = fCounterCollection->GetKeyWords("event").Tokenize(",");
 173
 174   Bool_t doPS = (events->FindObject("PSALL") != 0x0);
 175
 176   TObjString* srun;
 177   TObjString* strigger;
 178
 179   ULong64_t localNmb(0);
 180   ULong64_t localNmsl(0);
 181   ULong64_t localNmul(0);
 182
 183   if ( totalNmb) *totalNmb = 0;
 184   if ( totalNmsl) *totalNmsl = 0;
 185   if ( totalNmul ) *totalNmul = 0;
 186
 187   while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
 188   {
 189     std::cout << Form("RUN %09d ",srun->String().Atoi());
 190
 191     TString details;
 192     ULong64_t nmb(0);
 193     ULong64_t nmsl(0);
 194     ULong64_t nmul(0);
 195
 196     nextTrigger.Reset();
 197
 198     Int_t nofPS(0);
 199
 200     while ( ( strigger = static_cast<TObjString*>(nextTrigger()) ) )
 201     {
 202
 203       if ( !Config()->GetList(AliAnalysisMuMuConfig::kMinbiasTriggerList,IsSimulation()).Contains(strigger->String().Data()) &&
 204            !Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,IsSimulation()).Contains(strigger->String().Data()) &&
 205            !Config()->GetList(AliAnalysisMuMuConfig::kMuonTriggerList,IsSimulation()).Contains(strigger->String().Data()) ) continue;
 206
 207       ULong64_t n = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
 208                                                     strigger->String().Data(),"ALL",srun->String().Atoi())));
 209
 210       details += TString::Format("\n%50s %10lld",strigger->String().Data(),n);
 211
 212
 213       ULong64_t nps = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
 214                                                       strigger->String().Data(),"PSALL",srun->String().Atoi())));
 215
 216       if ( doPS )
 217       {
 218         details += TString::Format(" PS %5.1f %%",nps*100.0/n);
 219       }
 220
 221       if (nps)
 222       {
 223         ++nofPS;
 224       }
 225
 226       if ( Config()->GetList(AliAnalysisMuMuConfig::kMinbiasTriggerList,IsSimulation()).Contains(strigger->String()) )
 227       {
 228         nmb += n;
 229         if ( totalNmb) (*totalNmb) += n;
 230         localNmb += n;
 231       }
 232       else if ( Config()->GetList(AliAnalysisMuMuConfig::kMuonTriggerList,IsSimulation()).Contains(strigger->String()) )
 233       {
 234         nmsl += n;
 235         if ( totalNmsl) (*totalNmsl) += n;
 236         localNmsl += n;
 237       }
 238       else if ( Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,IsSimulation()).Contains(strigger->String()) )
 239       {
 240         nmul += n;
 241         if ( totalNmul ) (*totalNmul) += n;
 242         localNmul += n;
 243       }
 244     }
 245
 246     std::cout << Form("MB %10lld MSL %10lld MUL %10lld %s",
 247                  nmb,nmsl,nmul,(nofPS == 0 ? "(NO PS AVAIL)": ""));
 248
 249     if ( detailTriggers )
 250     {
 251       std::cout << details.Data();
 252     }
 253     std::cout << std::endl;
 254   }
 255
 256   if ( !totalNmul && !totalNmsl && !totalNmb )
 257   {
 258     std::cout << std::endl << Form("%13s MB %10lld MSL %10lld MUL %10lld ","TOTAL",
 259                                    localNmb,localNmsl,localNmul) << std::endl;
 260   }
 261
 262   delete runs;
 263   delete triggers;
 264   delete events;
 265 }
 266
 267
 268
 269 //_____________________________________________________________________________
 270 void AliAnalysisMuMu::CleanAllSpectra()
 271 {
 272   /// Delete all the spectra we may have
 273
 274   OC()->RemoveByType("AliAnalysisMuMuSpectra");
 275   Update();
 276 }
 277
 278
 279 //_____________________________________________________________________________
 280 TObjArray* AliAnalysisMuMu::CompareJpsiPerCMUUWithBackground(const char* jpsiresults,
 281                                                                    const char* backgroundresults)
 282 {
 283   TFile* fjpsi = FileOpen(jpsiresults);
 284   TFile* fbck = FileOpen(backgroundresults);
 285
 286   if (!fjpsi || !fbck) return 0x0;
 287
 288   TGraph* gjpsi = static_cast<TGraph*>(fjpsi->Get("jpsipercmuu"));
 289
 290   std::vector<std::string> checks;
 291
 292   checks.push_back("muminus-CMUU7-B-NOPF-ALLNOTRD");
 293   checks.push_back("muplus-CMUU7-B-NOPF-ALLNOTRD");
 294   checks.push_back("muminus-CMUSH7-B-NOPF-MUON");
 295   checks.push_back("muplus-CMUSH7-B-NOPF-MUON");
 296
 297   if (!gjpsi) return 0x0;
 298
 299   TObjArray* a = new TObjArray;
 300   a->SetOwner(kTRUE);
 301
 302   for ( std::vector<std::string>::size_type j = 0; j < checks.size(); ++j )
 303   {
 304
 305     TGraph* gback = static_cast<TGraph*>(fbck->Get(checks[j].c_str()));
 306
 307     if (!gback) continue;
 308
 309     if ( gjpsi->GetN() != gback->GetN() )
 310     {
 311       AliErrorClass("graphs have different number of points !");
 312       continue;
 313     }
 314
 315     TGraphErrors* g = new TGraphErrors(gjpsi->GetN());
 316
 317     for ( int i = 0; i < gjpsi->GetN(); ++i )
 318     {
 319       double r1,r2,y1,y2;
 320
 321       gjpsi->GetPoint(i,r1,y1);
 322       gback->GetPoint(i,r2,y2);
 323
 324       if ( r1 != r2 )
 325       {
 326         AliWarningClass(Form("run[%d]=%d vs %d",i,(int)r1,(int)r2));
 327         continue;
 328       }
 329
 330       g->SetPoint(i,y2,y1);
 331       //    g->SetPointError(i,gjpsi->GetErrorY(i),gback->GetErrorY(i));
 332     }
 333
 334     g->SetMarkerStyle(25+j);
 335     g->SetMarkerSize(1.2);
 336     if (j==0)
 337     {
 338       g->Draw("ap");
 339     }
 340     else
 341     {
 342       g->Draw("p");
 343     }
 344     g->SetLineColor(j+1);
 345     g->SetMarkerColor(j+1);
 346     g->SetName(checks[j].c_str());
 347     a->AddLast(g);
 348   }
 349
 350   return a;
 351 }
 352
 353 //_____________________________________________________________________________
 354 TGraph* AliAnalysisMuMu::CompareJpsiPerCMUUWithSimu(const char* realjpsiresults,
 355                                                              const char* simjpsiresults)
 356 {
 357   TFile* freal = FileOpen(realjpsiresults);
 358   TFile* fsim = FileOpen(simjpsiresults);
 359
 360   if (!freal || !fsim) return 0x0;
 361
 362   TGraph* greal = static_cast<TGraph*>(freal->Get("jpsipercmuu"));
 363   TGraph* gsim = static_cast<TGraph*>(fsim->Get("jpsipercmuu"));
 364
 365   TObjArray* a = new TObjArray;
 366   a->SetOwner(kTRUE);
 367
 368   if ( greal->GetN() != gsim->GetN() )
 369   {
 370     AliErrorClass("graphs have different number of points !");
 371     return 0x0;
 372   }
 373
 374   TGraphErrors* g = new TGraphErrors(greal->GetN());
 375   TGraphErrors* gratio = new TGraphErrors(greal->GetN());
 376
 377   for ( int i = 0; i < greal->GetN(); ++i )
 378   {
 379     double r1,r2,y1,y2;
 380
 381     greal->GetPoint(i,r1,y1);
 382     gsim->GetPoint(i,r2,y2);
 383
 384     if ( r1 != r2 )
 385     {
 386       AliWarningClass(Form("run[%d]=%d vs %d",i,(int)r1,(int)r2));
 387       continue;
 388     }
 389
 390     double ratio(0.0);
 391
 392     if ( TMath::Abs(y1)<1E-6 || TMath::Abs(y2)<1E-6)
 393     {
 394       g->SetPoint(i,0,0);
 395       g->SetPointError(i,0,0);
 396     }
 397     else
 398     {
 399       g->SetPoint(i,y2,y1);
 400       g->SetPointError(i,greal->GetErrorY(i),gsim ->GetErrorY(i));
 401       ratio = y2/y1;
 402     }
 403     gratio->SetPoint(i,r1,ratio);
 404   }
 405
 406   g->SetMarkerStyle(25);
 407   g->SetMarkerSize(1.2);
 408
 409   new TCanvas;
 410
 411   g->Draw("ap");
 412
 413   g->SetLineColor(1);
 414   g->SetMarkerColor(1);
 415   g->SetName("jpsipercmuurealvssim");
 416
 417   new TCanvas;
 418
 419   greal->Draw("alp");
 420   gsim->SetLineColor(4);
 421
 422   gsim->Draw("lp");
 423
 424   new TCanvas;
 425   gratio->Draw("alp");
 426
 427   return g;
 428 }
 429
 430 //_____________________________________________________________________________
 431 AliAnalysisMuMuConfig* AliAnalysisMuMu::Config()
 432 {
 433   /// Return the configuration
 434   return fConfig;
 435 }
 436
 437 //_____________________________________________________________________________
 438 void AliAnalysisMuMu::DrawMinv(const char* type,
 439                                const char* particle,
 440                                const char* trigger,
 441                                const char* eventType,
 442                                const char* pairCut,
 443                                const char* centrality,
 444                                const char* subresultname,
 445                                const char* flavour) const
 446 {
 447   /// Draw minv spectra for binning of given type
 448
 449   if (!OC() || !BIN()) return;
 450
 451   TObjArray* bins = BIN()->CreateBinObjArray(particle,type,flavour);
 452   if (!bins)
 453   {
 454     AliError(Form("Could not get %s bins",type));
 455     return;
 456   }
 457
 458   Double_t xmin(-1);
 459   Double_t xmax(-1);
 460
 461   TString sparticle(particle);
 462   if ( sparticle=="PSI" )
 463   {
 464     xmin = 2;
 465     xmax = 6;
 466   }
 467
 468   Int_t nx(1);
 469   Int_t ny(1);
 470
 471   Int_t n = bins->GetEntries();
 472
 473   if ( n == 2 )
 474   {
 475     nx = 2;
 476   }
 477   else if ( n > 2 )
 478   {
 479     ny = TMath::Nint(TMath::Sqrt(n));
 480     nx = n/ny;
 481   }
 482
 483   TString stype(type);
 484   stype.ToUpper();
 485
 486   TString spectraName(Form("/%s/%s/%s/%s/%s-%s",eventType,trigger,centrality,pairCut,particle,stype.Data()));
 487
 488   if ( strlen(flavour))
 489   {
 490     spectraName += "-";
 491     spectraName += flavour;
 492   }
 493
 494   AliAnalysisMuMuSpectra* spectra = static_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(spectraName.Data()));
 495
 496   AliDebug(1,Form("spectraName=%s spectra=%p",spectraName.Data(),spectra));
 497
 498   TObjArray* spectraBins(0x0);
 499   if ( spectra )
 500   {
 501     spectraBins = spectra->BinContentArray();
 502   }
 503
 504   TCanvas* c = new TCanvas;
 505   c->Divide(nx,ny);
 506   c->Draw();
 507   gStyle->SetOptFit(1112);
 508
 509   c->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", "AliAnalysisMuMu",
 510               (void*)this, "ExecuteCanvasEvent(Int_t,Int_t,Int_t,TObject*)");
 511
 512
 513   TIter next(bins);
 514   AliAnalysisMuMuBinning::Range* r;
 515   Int_t ci(0);
 516
 517   while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(next())) )
 518   {
 519     TString name(Form("/%s/%s/%s/%s/MinvUS%s",eventType,trigger,centrality,pairCut,r->AsString().Data()));
 520
 521     AliDebug(1,name.Data());
 522
 523     AliAnalysisMuMuJpsiResult* spectraBin(0x0);
 524
 525     if ( spectraBins )
 526     {
 527       AliAnalysisMuMuResult* sr = static_cast<AliAnalysisMuMuResult*>(spectraBins->At(ci));
 528
 529       spectraBin = static_cast<AliAnalysisMuMuJpsiResult*>(sr->SubResult(subresultname));
 530
 531       AliDebug(1,Form("spectraBin(%s)=%p",subresultname,spectraBin));
 532     }
 533
 534     TH1* h = OC()->Histo(name.Data());
 535
 536     if ( spectraBin )
 537     {
 538       h = spectraBin->Histo();
 539     }
 540
 541     if (h)
 542     {
 543       ++ci;
 544       c->cd(ci);
 545       gPad->SetLogy();
 546       if (xmin>0)
 547       {
 548         h->GetXaxis()->SetRangeUser(xmin,xmax);
 549       }
 550       h->Draw("histes");
 551
 552       TObject* f1 = h->GetListOfFunctions()->FindObject("fitTotal");
 553       if (f1)
 554       {
 555         f1->Draw("same");
 556       }
 557
 558       gPad->Modified();
 559       gPad->Update();
 560
 561       TObject* stats = h->FindObject("stats");
 562       if (stats)
 563       {
 564         stats->Draw("same");
 565       }
 566     }
 567   }
 568
 569   delete bins;
 570 }
 571
 572 //_____________________________________________________________________________
 573 void AliAnalysisMuMu::DrawMinv(const char* type, const char* particle, const char* flavour, const char* subresultname) const
 574 {
 575   /// Draw minv spectra for binning of given type
 576
 577 //  AliWarning("Reimplement me!");
 578
 579   if (!fConfig)
 580   {
 581     AliError("No configuration available yet. Don't know what to draw");
 582     return;
 583   }
 584
 585   DrawMinv(type,particle,
 586            First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,IsSimulation())),
 587            First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,IsSimulation())),
 588            First(Config()->GetList(AliAnalysisMuMuConfig::kPairSelectionList,IsSimulation())),
 589            First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,IsSimulation())),
 590            subresultname,
 591            flavour);
 592 //           First(Config()->kEventSelectionList()).Data(),
 593 //           First(Config()->kPairSelectionList()).Data(),
 594 //           First(Config()->kCentralitySelectionList()).Data(),
 595 //           subresultname,
 596 //           flavour);
 597 }
 598
 599 //___________________________________________________________________
 600 void AliAnalysisMuMu::ExecuteCanvasEvent(Int_t event, Int_t /*px*/, Int_t /*py*/, TObject *sel)
 601 {
 602   // Actions in reponse to mouse button events.
 603
 604   TCanvas* c = static_cast<TCanvas*>(gTQSender);
 605   TPad* pad = static_cast<TPad*>(c->GetSelectedPad());
 606   if (!pad) return;
 607
 608 //  if ((event == kButton1Down) ||
 609   if (event == kButton1Double)
 610   {
 611
 612 //    Float_t x = pad->AbsPixeltoX(px);
 613 //    Float_t y = pad->AbsPixeltoY(py);
 614 //    x = pad->PadtoX(x);
 615 //    y = pad->PadtoY(y);
 616
 617 //    std::cout << "event=" << event << " px=" << px << " py=" << py << " ";
 618
 619     if ( sel && sel->InheritsFrom("TH1") )
 620     {
 621       TCanvas* clocal = new TCanvas;
 622       clocal->SetLogy();
 623       clocal->Draw();
 624       sel->Draw();
 625     }
 626     else
 627     {
 628       TList* list = pad->GetListOfPrimitives();
 629       TIter next(list);
 630       TObject* h;
 631
 632       while ( ( h = next() ) )
 633       {
 634         if ( h->InheritsFrom("TH1") )
 635         {
 636           TCanvas* clocal = new TCanvas;
 637           clocal->SetLogy();
 638           clocal->Draw();
 639           h->Draw();
 640           break;
 641         }
 642       }
 643
 644     }
 645
 646 //      std::cout  << std::endl;
 647
 648       pad->Modified();
 649   }
 650
 651 }
 652
 653 //_____________________________________________________________________________
 654 TString
 655 AliAnalysisMuMu::ExpandPathName(const char* file)
 656 {
 657   // An expand method that lives alien URL as they are
 658   TString sfile;
 659
 660   if ( !sfile.BeginsWith("alien://") )
 661   {
 662     return gSystem->ExpandPathName(file);
 663   }
 664   else
 665   {
 666     if (!gGrid) TGrid::Connect("alien://");
 667     if (!gGrid) return "";
 668   }
 669
 670   return file;
 671 }
 672
 673 //_____________________________________________________________________________
 674 void AliAnalysisMuMu::TwikiOutputFnorm(const char* series) const
 675 {
 676   // make a twiki-compatible output of the Fnorm factor(s)
 677   TObjArray* what = TString(series).Tokenize(",");
 678   TObjString* s;
 679   TObjArray graphs;
 680   TIter next(what);
 681
 682   std::cout << "| *Run* |";
 683   while ( ( s = static_cast<TObjString*>(next())) )
 684   {
 685     TGraph* g = static_cast<TGraph*>(OC()->GetObject(Form("/FNORM/GRAPHS/%s",s->String().Data())));
 686     if (!g)
 687     {
 688       AliError(Form("Could not find graph for %s",s->String().Data()));
 689       continue;
 690     }
 691     std::cout << " *" << s->String().Data();
 692     if ( s->String().BeginsWith("RelDif") ) std::cout << " %";
 693     std::cout << "*|";
 694     graphs.Add(g);
 695   }
 696
 697   std::cout << std::endl;
 698
 699   TGraphErrors* g0 = static_cast<TGraphErrors*>(graphs.First());
 700   if (!g0) return;
 701
 702   for ( Int_t i = 0; i < g0->GetN(); ++i )
 703   {
 704     TString msg;
 705
 706     msg.Form("|%6d|",TMath::Nint(g0->GetX()[i]));
 707
 708     for ( Int_t j = 0; j < graphs.GetEntries(); ++j )
 709     {
 710       TGraphErrors* g = static_cast<TGraphErrors*>(graphs.At(j));
 711
 712       msg += TString::Format(" %6.2f +- %6.2f |",g->GetY()[i],g->GetEY()[i]);
 713     }
 714
 715     std::cout << msg.Data() << std::endl;
 716   }
 717
 718   next.Reset();
 719
 720   std::cout << "|*Weigthed mean (*)*|";
 721
 722   AliAnalysisMuMuResult* r = static_cast<AliAnalysisMuMuResult*>(OC()->GetObject("/FNORM/RESULTS/Fnorm"));
 723
 724   if (!r)
 725   {
 726     AliError("Could not find Fnorm result !");
 727     return;
 728   }
 729
 730
 731   while ( ( s = static_cast<TObjString*>(next())) )
 732   {
 733     TString var("Fnorm");
 734     TString unit;
 735
 736     if ( s->String().BeginsWith("Fnorm") )
 737     {
 738       r = static_cast<AliAnalysisMuMuResult*>(OC()->GetObject("/FNORM/RESULTS/Fnorm"));
 739     }
 740     else if ( s->String().BeginsWith("RelDif") )
 741     {
 742       r = static_cast<AliAnalysisMuMuResult*>(OC()->GetObject("/FNORM/RESULTS/RelDif"));
 743       unit = "%";
 744     }
 745
 746     r->Exclude("*");
 747     r->Include(s->String().Data());
 748
 749     std::cout << Form(" * %5.2f +- %5.2f %s * |",
 750                       r->GetValue(var.Data()),
 751                       r->GetErrorStat(var.Data()),
 752                       unit.Data());
 753   }
 754
 755   next.Reset();
 756
 757   std::cout << std::endl;
 758
 759   std::cout << "|*RMS*|";
 760
 761   while ( ( s = static_cast<TObjString*>(next())) )
 762   {
 763     TString var("Fnorm");
 764
 765     if ( s->String().BeginsWith("Fnorm") )
 766     {
 767       r = static_cast<AliAnalysisMuMuResult*>(OC()->GetObject("/FNORM/RESULTS/Fnorm"));
 768     }
 769     else if ( s->String().BeginsWith("RelDif") )
 770     {
 771       r = static_cast<AliAnalysisMuMuResult*>(OC()->GetObject("/FNORM/RESULTS/RelDif"));
 772     }
 773
 774     r->Exclude("*");
 775     r->Include(s->String().Data());
 776
 777     Double_t d = 100.0*r->GetRMS(var.Data())/r->GetValue(var.Data());
 778
 779     std::cout << Form(" * %5.2f (%5.2f %%) * |",
 780                       r->GetRMS(var.Data()),d);
 781   }
 782
 783   std::cout << std::endl;
 784   std::cout << "(*) weight is the number of CMUL7-B-NOPF-MUON triggers (physics-selected and pile-up corrected) in each run" << std::endl;
 785
 786   delete what;
 787 }
 788
 789 //_____________________________________________________________________________
 790 TFile*
 791 AliAnalysisMuMu::FileOpen(const char* file)
 792 {
 793   // Open a file after expansion of its name
 794
 795   return TFile::Open(ExpandPathName(file).Data());
 796 }
 797
 798 //_____________________________________________________________________________
 799 TString AliAnalysisMuMu::First(const TString& list) const
 800 {
 801   TObjArray* a = list.Tokenize(",");
 802   if ( a->GetLast() < 0 ) return "";
 803
 804   TString rv = static_cast<TObjString*>(a->First())->String();
 805
 806   delete a;
 807
 808   return rv;
 809 }
 810
 811 //_____________________________________________________________________________
 812 AliAnalysisMuMuSpectra*
 813 AliAnalysisMuMu::FitParticle(const char* particle,
 814                              const char* trigger,
 815                              const char* eventType,
 816                              const char* pairCut,
 817                              const char* centrality,
 818                              const AliAnalysisMuMuBinning& binning,
 819                              const char* spectraType,
 820                              Bool_t corrected)
 821 {
 822   // Fit the minv/mpt spectra to find the given particle
 823   // Returns an array of AliAnalysisMuMuResult objects
 824
 825   TProfile::Approximate(); //To avoid bins with error=0 due to low statstics
 826
 827   static int n(0);
 828
 829   TObjArray* bins = binning.CreateBinObjArray(particle);
 830   if (!bins)
 831   {
 832     AliError(Form("Did not get any bin for particle %s",particle));
 833     return 0x0;
 834   }
 835
 836   TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
 837   if ( !triggers->FindObject(trigger) )
 838   {
 839     AliError(Form("Did not find trigger %s",trigger));
 840     delete bins;
 841     delete triggers;
 842     return 0x0;
 843   }
 844   delete triggers;
 845
 846   TObjArray* events = fCounterCollection->GetKeyWords("event").Tokenize(",");
 847   if ( !events->FindObject(eventType) )
 848   {
 849     AliError(Form("Did not find eventType %s",eventType));
 850     delete bins;
 851     delete events;
 852     return 0x0;
 853   }
 854   delete events;
 855
 856   Int_t ntrigger = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s",trigger,eventType)));
 857
 858   if  (ntrigger<=0)
 859   {
 860     AliError(Form("No trigger for trigger:%s/event:%s",trigger,eventType));
 861     delete bins;
 862     return 0x0;
 863   }
 864
 865   TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
 866   Int_t nruns = runs->GetEntries();
 867   delete runs;
 868
 869
 870   TString id(Form("/%s/%s/%s/%s",eventType,trigger,centrality,pairCut));
 871
 872 //  binning.Print();
 873
 874   AliAnalysisMuMuSpectra* spectra(0x0);
 875
 876   AliAnalysisMuMuBinning::Range* bin;
 877   TIter next(bins);
 878
 879   TObjArray* fitTypeArray = Config()->GetList(AliAnalysisMuMuConfig::kFitTypeList,IsSimulation()).Tokenize(",");
 880   TIter nextFitType(fitTypeArray);
 881   TObjString* fitType;
 882   TString flavour;
 883   TString sSpectraType(spectraType);
 884
 885   TString spectraName(binning.GetName());
 886   if ( flavour.Length() > 0 )
 887   {
 888     spectraName += "-";
 889     spectraName += flavour;
 890   }
 891   if ( corrected )
 892   {
 893     spectraName += "-";
 894     spectraName += "AccEffCorr";
 895   }
 896
 897 //  Int_t binN(0);
 898   while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(next())) )
 899   {
 900     TString hname;
 901     if (!sSpectraType.CompareTo("minv"))
 902     {
 903       hname = corrected ? Form("MinvUS%s_AccEffCorr",bin->AsString().Data()) : Form("MinvUS%s",bin->AsString().Data());
 904     }
 905     else if (!sSpectraType.CompareTo("mpt"))
 906     {
 907       hname = corrected ? Form("MeanPtVsMinvUS%s_AccEffCorr",bin->AsString().Data()) : Form("MeanPtVsMinvUS%s",bin->AsString().Data());
 908     }
 909     else
 910     {
 911       AliError("Wrong spectra type choice: Posibilities are: 'minv' or 'mpt' ");
 912       return 0x0;
 913     }
 914
 915     TString isCorr(corrected ? " AccEffCorr " : " ");
 916     std::cout << "---------------------------------//---------------------------------" << std::endl;
 917     std::cout << "Fitting" << isCorr.Data() << sSpectraType.Data() << " spectra in " << id.Data() << std::endl;
 918
 919     TH1* histo = OC()->Histo(id.Data(),hname.Data());
 920
 921     if (!histo)
 922     {
 923 //      if (!fBinning && bin->IsIntegrated() )
 924 //      {
 925 //        // old file, we only had MinvUSPt
 926 //        hminv = fMergeableCollection->Histo(Form("/%s/%s/%s/%s",eventType,trigger,centrality,pairCut),"MinvUSPt:py");
 927 //      }
 928 //
 929 //      if (!hminv)
 930 //      {
 931         AliError(Form("Could not find histo %s",hname.Data()));
 932         continue;
 933 //      }
 934     }
 935
 936     histo = static_cast<TH1*>(histo->Clone(Form("%s%d",sSpectraType.Data(),n++)));
 937
 938     const char* particleTmp = IsSimulation() ? GetParticleName() : "JPsi"; //At some point particleTmp should become particle (but for now particle is always = "psi")
 939
 940     TString sparticleTmp(particleTmp);
 941
 942     AliAnalysisMuMuJpsiResult* r = new AliAnalysisMuMuJpsiResult(particleTmp,
 943                                                                  *histo, // Result for current bin
 944                                                                  trigger,
 945                                                                  eventType,
 946                                                                  pairCut,
 947                                                                  centrality,
 948                                                                  *bin);
 949
 950     r->SetNofTriggers(ntrigger);
 951     r->SetNofRuns(nruns);
 952
 953     nextFitType.Reset();
 954
 955     Int_t added(0);
 956
 957     while ( ( fitType = static_cast<TObjString*>(nextFitType())) )
 958     {
 959       // In this loop we create a Subresult for each fit inside the Result for current bin (AddFit will do)
 960
 961       TString sFitType(fitType->String());
 962
 963       if ( !sFitType.Contains(sSpectraType.Data()) ) continue;
 964
 965       AliDebug(1,Form("<<<<<< fitType=%s bin=%s",sFitType.Data(),bin->Flavour().Data()));
 966
 967       std::cout << "" << std::endl;
 968       std::cout << "---------------" << "Fit " << added + 1 << "------------------" << std::endl;
 969       std::cout << "Fitting " << hname.Data() << " with " << sFitType.Data() << std::endl;
 970       std::cout << "" << std::endl;
 971
 972       if ( sFitType.Contains("mctails",TString::kIgnoreCase) ) //FIXME: Find a univoque way to determine the correctly the fit type
 973       {
 974         TString sbin = bin->AsString();
 975         TString spectraMCName = spectraName;
 976         AliAnalysisMuMuBinning::Range* binMC = bin;
 977
 978         if ((sbin.Contains("MULT") || sbin.Contains("NCH") || sbin.Contains("DNCHDETA") || sbin.Contains("V0A") || sbin.Contains("V0ACENT") || sbin.Contains("NTRCORR")|| sbin.Contains("RELNTRCORR")) && !sbin.Contains("NTRCORRPT") && !sbin.Contains("NTRCORRY"))
 979         {
 980           //-------has to have a better way to do it
 981           AliAnalysisMuMuBinning* b = new AliAnalysisMuMuBinning;
 982           b->AddBin("psi","INTEGRATED");
 983
 984           binMC = static_cast<AliAnalysisMuMuBinning::Range*>(b->CreateBinObjArray()->At(0));
 985
 986           spectraMCName = b->GetName();
 987           delete b;
 988
 989 //          if ( flavour.Length() > 0 ) //Commented cause we set no flavour in "INTEGRATED" bin in the analysis task
 990 //          {
 991 //            spectraMCName += "-";
 992 //            spectraMCName += flavour;
 993 //          }
 994           if ( corrected )
 995           {
 996             spectraMCName += "-";
 997             spectraMCName += "AccEffCorr";
 998           }
 999           //-----------
1000
1001         }
1002 //        if( sbin.Contains("NTRCORRPT") || !sbin.Contains("NTRCORRY") )
1003 //        {
1004 //          spectraMCName.Remove(4,7);
1005 //
1006 //          if ( spectraMCName.Contains("PT") )
1007 //          {
1008 //            AliAnalysisMuMuBinning* b = SIM()->BIN();
1009 //
1010 //            TObjArray* binsPt = b->CreateBinObjArray(particle,"PT","");
1011 //
1012 //            Int_t nEntrBinMC = binsPt->GetEntries();
1013 //
1014 ////            binsPt->Print();
1015 //
1016 //            binMC = static_cast<AliAnalysisMuMuBinning::Range*>(binsPt->At(binN));
1017 //
1018 //            if ( binN == nEntrBinMC - 1 ) binN = -1;
1019 //
1020 //          }
1021 //
1022 //        }
1023
1024         //par = GetCB2Tails(*binInt,"MCTAILS",eventType,trigger,pairCut,corrected);Why I was taking the tails from the fitted data spectra?
1025         GetParametersFromMC(sFitType,Form("/%s/%s",centrality,pairCut),spectraMCName.Data(),binMC);
1026
1027         if (sFitType.Length()>0)
1028         {
1029           added += ( r->AddFit(sFitType.Data()) == kTRUE );
1030         }
1031       }
1032
1033       else if ( sFitType.Contains("mpt",TString::kIgnoreCase) && !sFitType.Contains("minv",TString::kIgnoreCase) )
1034       {
1035         std::cout << "++The Minv parameters will be taken from " << spectraName.Data() << std::endl;
1036         std::cout << "" << std::endl;
1037
1038         AliAnalysisMuMuSpectra* minvSpectra = dynamic_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(id.Data(),spectraName.Data()));
1039
1040         if ( !minvSpectra )
1041         {
1042           AliError(Form("Cannot fit mean pt: could not get the minv spectra for %s",id.Data()));
1043           continue;//return 0x0;
1044         }
1045
1046         AliAnalysisMuMuJpsiResult* minvResult = static_cast<AliAnalysisMuMuJpsiResult*>(minvSpectra->GetResultForBin(*bin));
1047
1048         if ( !minvResult )
1049         {
1050           AliError(Form("Cannot fit mean pt: could not get the minv result for bin %s in %s",bin->AsString().Data(),id.Data()));
1051           continue; //return 0x0;
1052         }
1053
1054         TObjArray* minvSubResults = minvResult->SubResults();
1055         TIter nextSubResult(minvSubResults);
1056         AliAnalysisMuMuJpsiResult* fitMinv;
1057         TString subResultName;
1058
1059         Int_t nSubFit(0);
1060         while ( ( fitMinv = static_cast<AliAnalysisMuMuJpsiResult*>(nextSubResult())) )
1061         {
1062           std::cout << "" << std::endl;
1063           std::cout <<  "      /-- SubFit " << nSubFit + 1 << " --/ " << std::endl;
1064           std::cout << "" << std::endl;
1065
1066           TString fitMinvName(fitMinv->GetName());
1067           fitMinvName.Remove(fitMinvName.First("_"),fitMinvName.Sizeof()-fitMinvName.First("_"));
1068
1069 //          std::cout << fitMinvName.Data() << " ; " << fitMinv->GetName() << std::endl;
1070 //          std::cout << "" << std::endl;
1071
1072           if ( !sFitType.Contains(fitMinvName) ) continue; //FIXME: Ambiguous, i.e. NA60NEWPOL2EXP & NA60NEWPOL2 (now its ok cause only VWG and POL2EXP are used, but care)
1073
1074           TString sMinvFitType(sFitType);
1075
1076           GetParametersFromResult(sMinvFitType,fitMinv);//FIXME: Think about if this is necessary
1077
1078           added += ( r->AddFit(sMinvFitType.Data()) == kTRUE );
1079
1080           nSubFit++;
1081         }
1082       }
1083
1084       else if ( sFitType.Contains("minv&mpt",TString::kIgnoreCase) )
1085       {
1086         AliWarning("Implement here the method to do the combined minv mpt fits");
1087         //FIXME: Shall we use the fitType or spectraType to choose to perform combined fits? Cause we have to check what kind of object is returned by the combined fit in order to decide if we put it in a different spectra(spectraType would be the flag,and here we should update the spectraName) or as a subresult(fitType in this case)
1088       }
1089
1090       else
1091       {
1092         if ( sFitType.Contains("PSICB2",TString::kIgnoreCase) || sFitType.Contains("PSINA60NEW",TString::kIgnoreCase)) std::cout << "+Free tails fit... " << std::endl;
1093         else if ( sFitType.Contains("PSICOUNT",TString::kIgnoreCase) )  std::cout << Form("+Just counting %s...",GetParticleName()) << std::endl;
1094         else std::cout << "+Using predefined tails... " << std::endl;
1095
1096         if ( sFitType.Contains("minvJPsi") && !sparticleTmp.Contains("JPsi") )
1097         {
1098           std::cout << "This fitting funtion is set to fit JPsi: Skipping fit..." << std::endl;
1099           continue;
1100         }
1101         if ( sFitType.Contains("minvPsiP") && !sparticleTmp.Contains("PsiP") )
1102         {
1103           std::cout << "This fitting funtion is set to fit PsiP: Skipping fit..." << std::endl;
1104           continue;
1105         }
1106
1107         added += ( r->AddFit(sFitType.Data()) == kTRUE );
1108       }
1109
1110       std::cout << "-------------------------------------" << std::endl;
1111       std::cout << "" << std::endl;
1112     }
1113
1114     if ( !added ) continue;
1115
1116 //    TObjArray* a = r->SubResults(); // TEST
1117 //    a->Print();
1118
1119     flavour = bin->Flavour();
1120
1121     if (!spectra)
1122     {
1123       TString spectraSaveName = spectraName;
1124       if ( !sSpectraType.CompareTo("mpt") )
1125       {
1126         spectraSaveName += "-";
1127         spectraSaveName += "MeanPtVsMinvUS";
1128       }
1129
1130       spectra = new AliAnalysisMuMuSpectra(spectraSaveName.Data());
1131     }
1132
1133     Bool_t adoptOk = spectra->AdoptResult(*bin,r); // We adopt the Result for current bin into the spectra
1134
1135     if ( adoptOk ) std::cout << "Result " << r->GetName() << " adopted in spectra " << spectra->GetName() << std::endl;
1136     else AliError(Form("Error adopting result %s in spectra %s",r->GetName(),spectra->GetName()));
1137
1138
1139     if ( IsSimulation() )
1140     {
1141       SetNofInputParticles(*r);
1142     }
1143
1144 //    binN++;
1145   } // loop on bins
1146
1147   delete fitTypeArray;
1148   delete bins;
1149
1150   return spectra;
1151 }
1152
1153
1154 ////_____________________________________________________________________________
1155 //AliAnalysisMuMuSpectra*
1156 //AliAnalysisMuMu::FitMeanPt(const char* particle,
1157 //                           const char* trigger,
1158 //                           const char* eventType,
1159 //                           const char* pairCut,
1160 //                           const char* centrality,
1161 //                           const AliAnalysisMuMuBinning& binning,const AliAnalysisMuMuSpectra& spectraMinv,Bool_t corrected)
1162 //{
1163 //  // Fit the dimuon mean pt spectra to find the particle mean pt
1164 //  // Returns an array of AliAnalysisMuMuResult objects
1165 //
1166 //  TObjArray* bins = binning.CreateBinObjArray(particle);
1167 //  if (!bins)
1168 //  {
1169 //    AliError(Form("Did not get any bin for particle %s",particle));
1170 //    return 0x0;
1171 //  }
1172 //
1173 //  //  AliAnalysisMuMuBinning* binningMinv(0x0);
1174 //  //  TObjArray* binsMinv(0x0);
1175 //  //  AliAnalysisMuMuBinning::Range* binInt(0x0);
1176 //  //  AliAnalysisMuMuBinning::Range* testbin = static_cast<AliAnalysisMuMuBinning::Range*>(bins->At(0));
1177 //  //  TString sbins = testbin->AsString();
1178 //
1179 //  //  if ( sbins.Contains("MULT")) //For the multiplicity bins we will use the tails from the integrated Minv spectra
1180 //  //  {
1181 //  //    binningMinv = spectraMinv.Binning();
1182 //  //    binsMinv = binningMinv->CreateBinObjArray(particle);
1183 //  //    binInt = static_cast<AliAnalysisMuMuBinning::Range*>(binsMinv->At(0));
1184 //  //  }
1185 //
1186 //
1187 //  TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
1188 //  if ( !triggers->FindObject(trigger) )
1189 //  {
1190 //    AliDebug(1,Form("Did not find trigger %s",trigger));
1191 //    delete bins;
1192 //    delete triggers;
1193 //    return 0x0;
1194 //  }
1195 //  delete triggers;
1196 //
1197 //  TObjArray* events = fCounterCollection->GetKeyWords("event").Tokenize(",");
1198 //  if ( !events->FindObject(eventType) )
1199 //  {
1200 //    AliError(Form("Did not find eventType %s",eventType));
1201 //    delete bins;
1202 //    delete events;
1203 //    return 0x0;
1204 //  }
1205 //  delete events;
1206 //
1207 //  Int_t ntrigger = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s",trigger,eventType)));
1208 //
1209 //  if  (ntrigger<=0)
1210 //  {
1211 //    AliError(Form("No trigger for trigger:%s/event:%s",trigger,eventType));
1212 //    delete bins;
1213 //    return 0x0;
1214 //  }
1215 //
1216 //  AliAnalysisMuMuSpectra* sMinv = static_cast<AliAnalysisMuMuSpectra*>(spectraMinv.Clone());
1217 //  if(!sMinv)
1218 //  {
1219 //    AliError("Did not find inv mass spectra");
1220 //  }
1221 //
1222 //
1223 //  //  binning.Print();
1224 //
1225 //  AliAnalysisMuMuSpectra* spectra(0x0);
1226 //
1227 //  AliAnalysisMuMuBinning::Range* bin;
1228 //  //  AliAnalysisMuMuBinning::Range* binMinv;
1229 //  TIter next(bins);
1230 //
1231 //  //  TObjArray* fitTypeArray = fFitTypeList.Tokenize(",");
1232 //  //  TIter nextFitType(fitTypeArray);
1233 //  //  TObjString* fitType;
1234 //  TString flavour;
1235 //
1236 //
1237 //  while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(next())) )
1238 //  {
1239 //    std::cout << "---------------------------------" << std::endl;
1240 //    std::cout << "Fitting mean Pt in " << bin->AsString().Data() << " " << "for" << " " << eventType << "/" << trigger << "/" << centrality << "/" << pairCut << std::endl;
1241 //    std::cout << "---------------------------------" << std::endl;
1242 //
1243 //
1244 //    TString pname;
1245 //    if ( corrected ) pname = Form("MeanPtVsMinvUS%s_Corr",bin->AsString().Data());
1246 //    else pname = Form("MeanPtVsMinvUS%s",bin->AsString().Data());
1247 //
1248 //    TProfile* hmPt = static_cast<TProfile*>(fMergeableCollection->GetObject(Form("/%s/%s/%s/%s",eventType,trigger,centrality,pairCut),pname.Data()));
1249 //
1250 //    if (!hmPt)
1251 //    {
1252 //      //      if (!fBinning && bin->IsNullObject() )
1253 //      //      {
1254 //      //        // old file, we only had MinvUSPt
1255 //      //        hminv = fMergeableCollection->Histo(Form("/%s/%s/%s/%s",eventType,trigger,centrality,pairCut),"MinvUSPt:py");
1256 //      //      }
1257 //
1258 //      //      if (!hmPt)
1259 //      //      {
1260 //      AliDebug(1,Form("Could not find histo profile %s",pname.Data()));
1261 //      continue;
1262 //      //      }
1263 //    }
1264 //    hmPt->Approximate();
1265 //
1266 //
1267 //    //    hmPt = static_cast<TH1*>(hmPt->Clone(Form("mPtv%d",n++))); //Reuse this
1268 //
1269 //    //    if ( sbins.Contains("MULT") )
1270 //    //    {
1271 //    //      binMinv = binInt;
1272 //    //    }
1273 //    //    else
1274 //    //    {
1275 //    //      binMinv = bin;
1276 //    //    }
1277 //
1278 //    AliAnalysisMuMuJpsiResult* rMinv = static_cast<AliAnalysisMuMuJpsiResult*>(sMinv->GetResultForBin(*bin/*Minv*/));
1279 //
1280 //    if( !rMinv )
1281 //    {
1282 //      AliError(Form("Did not find inv mass result inside spectra for bin %s",bin/*Minv*/->Quantity().Data()));
1283 //    }
1284 //
1285 //    TObjArray* fits = rMinv->SubResults();
1286 //    AliAnalysisMuMuResult* fitMinv;
1287 //    TIter nextFit(fits);
1288 //    TString fitName;
1289 //
1290 //    AliAnalysisMuMuJpsiResult* r = new AliAnalysisMuMuJpsiResult(*hmPt,
1291 //                                                         trigger,
1292 //                                                         eventType,
1293 //                                                         pairCut,
1294 //                                                         centrality,
1295 //                                                         *bin);
1296 //
1297 //    r->SetNofTriggers(ntrigger);
1298 //
1299 //    nextFit.Reset();
1300 //
1301 //    //    new TCanvas;
1302 //    Int_t i = 1;
1303 //
1304 //    while ( ( fitMinv = static_cast<AliAnalysisMuMuJpsiResult*>(nextFit())) )
1305 //    {
1306 //      std::cout << "" << std::endl;
1307 //      std::cout << "---------------" << "Fit " << i << "------------------" << std::endl;
1308 //      i++;
1309 //
1310 //      fitName = fitMinv->Alias();
1311 //      std::cout << "" << std::endl;
1312 //      std::cout << "Getting signal parameters from " << eventType << "/" << trigger << "/" << centrality << "/" << pairCut << "/" << spectraMinv.GetName()  << std::endl;
1313 //      std::cout << "" << std::endl;
1314 //
1315 //      if(fitName.BeginsWith("JPSI2CB2VWG") || fitName.BeginsWith("MCTAILS"))
1316 //      {
1317 //        Double_t par[14] = {fitMinv->GetValue("kVWG"),fitMinv->GetValue("mVWG"),fitMinv->GetValue("sVWG1"),
1318 //          fitMinv->GetValue("sVWG2"),fitMinv->GetValue("kPsi"),fitMinv->GetValue("MeanJpsi"),fitMinv->GetValue("SigmaJpsi"),
1319 //          fitMinv->GetValue("alPsi"),fitMinv->GetValue("nlPsi"),fitMinv->GetValue("auPsi"),fitMinv->GetValue("nuPsi"),fitMinv->GetValue("kPsi'"),fitMinv->GetValue("NofJpsi"),fitMinv->GetErrorStat("NofJpsi")}; //Create a method to get the parameters //The last 2 parameters are not used in the fit
1320 //
1321 //        if(fitName.BeginsWith("JPSI2CB2VWG"))
1322 //        {
1323 //          std::cout << "PRE-SET TAILS" << std::endl;
1324 //          std::cout << "" << std::endl;
1325 //          r->AddMeanPtFit("JPSI2CB2VWG:2",&par[0]); // FIXME: :2 should be in the default fit types but for meanpt fit naming is not there, change it
1326 //        }
1327 //        else
1328 //        {
1329 //          std::cout << "MC TAILS" << std::endl;
1330 //          std::cout << "" << std::endl;
1331 //          r->AddMeanPtFit("MCTAILS-JPSI2CB2VWG:2",&par[0]);
1332 //        }
1333 //
1334 //
1335 //      }
1336 //
1337 //      //Make a part for NA48 and the rest of fitting functions.
1338 //
1339 //    }
1340 //
1341 //    flavour = bin->Flavour();
1342 //
1343 //    if (!spectra)
1344 //    {
1345 //      TString spectraName(Form("MeanPtVsMinvUS-%s",binning.GetName()));
1346 //      if ( flavour.Length() > 0 )
1347 //      {
1348 //        spectraName += "-";
1349 //        spectraName += flavour;
1350 //      }
1351 //      if ( corrected )
1352 //      {
1353 //        spectraName += "-";
1354 //        spectraName += "Corr";
1355 //      }
1356 //
1357 //      spectra = new AliAnalysisMuMuSpectra(spectraName.Data());
1358 //    }
1359 //
1360 //    spectra->AdoptResult(*bin,r);
1361 //
1362 //    //    if ( IsSimulation() )
1363 //    //    {
1364 //    //      SetNofInputParticles(*r);
1365 //    //    }
1366 //
1367 //
1368 //  } // loop on bins
1369 //
1370 //  //  delete fitTypeArray;
1371 //  delete sMinv;
1372 //  delete bins;
1373 //
1374 //  return spectra;
1375 //}
1376
1377 //_____________________________________________________________________________
1378 void
1379 AliAnalysisMuMu::GetParametersFromMC(TString& fitType, const char* pathCentrPairCut, const char* spectraName, AliAnalysisMuMuBinning::Range* bin) const
1380 {
1381   /// Get the parameters from the associated simulation. Is intended to be used in minv fits, where we need the tails from JPsi (and Psi')
1382   // We can choose between the 2 associated simulations (a JPsi and Psi' ones) by setting the sim variable to "sim1" (fAssociatedSimulation by default) or "sim2" (fAssociatedSimulation2)
1383
1384   if ( !SIM() && !SIM2() )
1385   {
1386     AliError("Cannot get MC tails without associated simulation(s) file(s) !");
1387     fitType = "";
1388     return;
1389   }
1390
1391   TString subResultName("");
1392   if ( fitType.Contains("NA60NEW",TString::kIgnoreCase) ) subResultName = "PSINA60NEW";//_1.5_4.2
1393   else if ( fitType.Contains("CB2",TString::kIgnoreCase) ) subResultName = "PSICB2";//_2.2_3.9
1394   else
1395   {
1396     AliError("I don't know from which MC subresult take the tails");
1397     return;
1398   }
1399
1400   TObjArray* simArr = new TObjArray;
1401   if ( SIM() ) simArr->Add(SIM());
1402   if ( SIM2() && fitType.Contains("INDEPTAILS",TString::kIgnoreCase) ) simArr->Add(SIM2()); // If we have set the key to get the JPsi ans PsiP tails
1403
1404   TIter nextSim(simArr);
1405   AliAnalysisMuMu* currentSIM;
1406
1407   TString spath(pathCentrPairCut);
1408
1409   spath.Prepend(Form("/%s",First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kTRUE)).Data()));//FIXME: Care with this when there is more than one selection in the list
1410   spath.Prepend(Form("/%s",First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kTRUE)).Data()));
1411
1412
1413   while ( (currentSIM = static_cast<AliAnalysisMuMu*>(nextSim())) )
1414   {
1415     AliAnalysisMuMuSpectra* minvMCSpectra = currentSIM->SPECTRA(Form("%s/%s",spath.Data(),spectraName));
1416     if (!minvMCSpectra)
1417     {
1418       AliError(Form("Could not find spectra %s/%s for associated simulation",spath.Data(),spectraName));
1419       currentSIM->OC()->Print("*:Ali*");
1420       fitType = "";
1421       continue;
1422     }
1423
1424     AliAnalysisMuMuJpsiResult* minvMCResult = static_cast<AliAnalysisMuMuJpsiResult*>(minvMCSpectra->GetResultForBin(*bin));
1425     if ( !minvMCResult )
1426     {
1427       AliError(Form("Cannot get MC tails cause the minv result for bin %s in %s/%s was not found",bin->AsString().Data(),spath.Data(),spectraName));
1428       fitType = "";
1429       continue;
1430     }
1431
1432     AliAnalysisMuMuJpsiResult* r = dynamic_cast<AliAnalysisMuMuJpsiResult*>(minvMCResult->SubResult(subResultName.Data())); //FIXME: Find an independet way of naming results
1433     if  ( r && subResultName.Contains("PSICB2") )
1434     {
1435       fitType += Form(":al%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("al%s",currentSIM->GetParticleName())));
1436       fitType += Form(":nl%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("nl%s",currentSIM->GetParticleName())));
1437       fitType += Form(":au%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("au%s",currentSIM->GetParticleName())));
1438       fitType += Form(":nu%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("nu%s",currentSIM->GetParticleName())));
1439
1440       std::cout << " Using MC " << currentSIM->GetParticleName() << " CB2 tails... " << std::endl;
1441       std::cout << std::endl;
1442     }
1443     else if ( r && subResultName.Contains("PSINA60NEW") )
1444     {
1445       fitType += Form(":p1L%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("p1L%s",currentSIM->GetParticleName())));
1446       fitType += Form(":p2L%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("p2L%s",currentSIM->GetParticleName())));
1447       fitType += Form(":p3L%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("p3L%s",currentSIM->GetParticleName())));
1448       fitType += Form(":p1R%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("p1R%s",currentSIM->GetParticleName())));
1449       fitType += Form(":p2R%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("p2R%s",currentSIM->GetParticleName())));
1450       fitType += Form(":p3R%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("p3R%s",currentSIM->GetParticleName())));
1451
1452       fitType += Form(":aL%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("aL%s",currentSIM->GetParticleName())));
1453       fitType += Form(":aR%s=%f",currentSIM->GetParticleName(),r->GetValue(Form("aR%s",currentSIM->GetParticleName())));
1454
1455       std::cout << " Using MC " << currentSIM->GetParticleName() << " NA60New tails... " << std::endl;
1456       std::cout << std::endl;
1457     }
1458     else
1459     {
1460       AliError(Form("Cannot get MC tails. MC Subresult %s not found",minvMCResult->GetName()));
1461       fitType = "";
1462     }
1463   }
1464   delete simArr;
1465 }
1466
1467 //_____________________________________________________________________________
1468 void
1469 AliAnalysisMuMu::GetParametersFromResult(TString& fitType, AliAnalysisMuMuJpsiResult* minvResult) const
1470 {
1471   // Gets the parameters from a result, is intended to be used for mean pt fits where we need the signal and backgroud parameters
1472
1473   AliWarning("Re-implement me !!!"); //FIXME: The parameters to get will depend on the fit function and also in this way is not suitable for other particles (ie Upsilon)(Find a way to get the particle(s) name)
1474
1475   TString msg("");
1476   if ( minvResult )
1477   {
1478     // We take the usual parameters (the ones from JPsi and the normalization of the Psi')
1479     fitType += Form(":kJPsi=%f",minvResult->GetValue("kJPsi")); //FIXME: Names are not correct
1480     fitType += Form(":mJPsi=%f",minvResult->GetValue("mJPsi"));
1481     fitType += Form(":sJPsi=%f",minvResult->GetValue("sJPsi"));
1482
1483     fitType += Form(":NofJPsi=%f",minvResult->GetValue("NofJPsi"));
1484     fitType += Form(":ErrStatNofJPsi=%f",minvResult->GetErrorStat("NofJPsi"));
1485
1486     fitType += Form(":kPsiP=%f",minvResult->GetValue("kPsiP"));
1487
1488     TString minvName(minvResult->GetName());
1489
1490     TString minvRangeParam = minvName;
1491     minvRangeParam.Remove(0,minvRangeParam.First("_") + 1);
1492     fitType += Form(":MinvRS=%s",minvRangeParam.Data());
1493
1494     fitType += Form(":FSigmaPsiP=%f",minvResult->GetValue("FSigmaPsiP"));
1495
1496     if ( fitType.Contains("CB2",TString::kIgnoreCase) )
1497     {
1498       fitType += Form(":alJPsi=%f",minvResult->GetValue("alJPsi"));
1499       fitType += Form(":nlJPsi=%f",minvResult->GetValue("nlJPsi"));
1500       fitType += Form(":auJPsi=%f",minvResult->GetValue("auJPsi"));
1501       fitType += Form(":nuJPsi=%f",minvResult->GetValue("nuJPsi"));
1502
1503       msg += "JPsi CB2 signal parameters";
1504       //    fitType += Form(":NofPsiP=%f",minvResult->GetValue("NofPsiP"));
1505       //    fitType += Form(":ErrStatNofPsiP=%f",minvResult->GetErrorStat("NofPsiP"));
1506
1507       if ( fitType.Contains("INDEPTAILS") )
1508       {
1509 //        minvName = minvResult->GetName();
1510         if ( minvName.Contains("INDEPTAILS") )
1511         {
1512           // In case we use independent parameters tails for JPsi and Psi' we take also the Psi' ones
1513           fitType += Form(":alPsiP=%f",minvResult->GetValue("alPsiP"));
1514           fitType += Form(":nlPsiP=%f",minvResult->GetValue("nlPsiP"));
1515           fitType += Form(":auPsiP=%f",minvResult->GetValue("auPsiP"));
1516           fitType += Form(":nuPsiP=%f",minvResult->GetValue("nuPsiP"));
1517           fitType += Form(":mPsiP=%f",minvResult->GetValue("mPsiP"));
1518           fitType += Form(":sPsiP=%f",minvResult->GetValue("sPsiP"));
1519
1520           msg += " + PsiP CB2 signal parameters";
1521         }
1522         else
1523         {
1524           AliError(Form("Cannot get PsiP tails from result. Result %s does not contain PsiP tails info => Fit will fail ",minvResult->GetName()));
1525           fitType = "";
1526           return;
1527         }
1528       }
1529     }
1530     else if ( fitType.Contains("NA60NEW",TString::kIgnoreCase) )
1531     {
1532       fitType += Form(":p1LJPsi=%f",minvResult->GetValue("p1LJPsi"));
1533       fitType += Form(":p2LJPsi=%f",minvResult->GetValue("p2LJPsi"));
1534       fitType += Form(":p3LJPsi=%f",minvResult->GetValue("p3LJPsi"));
1535       fitType += Form(":p1RJPsi=%f",minvResult->GetValue("p1RJPsi"));
1536       fitType += Form(":p2RJPsi=%f",minvResult->GetValue("p2RJPsi"));
1537       fitType += Form(":p3RJPsi=%f",minvResult->GetValue("p3RJPsi"));
1538
1539       fitType += Form(":aLJPsi=%f",minvResult->GetValue("aLJPsi"));
1540       fitType += Form(":aRJPsi=%f",minvResult->GetValue("aRJPsi"));
1541
1542       msg += "JPsi NA60New signal parameters";
1543
1544       if ( fitType.Contains("INDEPTAILS") )
1545       {
1546 //        TString minvName(minvResult->GetName());
1547         if ( minvName.Contains("INDEPTAILS") )
1548         {
1549           // In case we use independent parameters tails for JPsi and Psi' we take also the Psi' ones
1550           fitType += Form(":p1LPsiP=%f",minvResult->GetValue("p1LPsiP"));
1551           fitType += Form(":p2LPsiP=%f",minvResult->GetValue("p2LPsiP"));
1552           fitType += Form(":p3LPsiP=%f",minvResult->GetValue("p3LPsiP"));
1553           fitType += Form(":p1RPsiP=%f",minvResult->GetValue("p1RPsiP"));
1554           fitType += Form(":p2RPsiP=%f",minvResult->GetValue("p2RPsiP"));
1555           fitType += Form(":p3RPsiP=%f",minvResult->GetValue("p3RPsiP"));
1556
1557           fitType += Form(":aLPsiP=%f",minvResult->GetValue("aLPsiP"));
1558           fitType += Form(":aRPsiP=%f",minvResult->GetValue("aRPsiP"));
1559
1560           msg += " + PsiP NA60New signal parameters";
1561
1562         }
1563         else
1564         {
1565           AliError(Form("Cannot get PsiP tails from result. Result %s does not contain PsiP tails info => Fit will fail ",minvResult->GetName()));
1566           fitType = "";
1567           return;
1568         }
1569       }
1570     }
1571     else
1572     {
1573       AliError(Form("Cannot get the parameters from %s",minvResult->GetName()));
1574       fitType = "";
1575       return;
1576     }
1577     // Now we take the background parameters
1578     if ( fitType.Contains("VWG_") || fitType.Contains("VWGINDEPTAILS") ) //FIXME: Check that cannot be misunderstood(like Exp x Pol2..). In fact it can be misunderstood since the meanpt function name has also the name of the function to fit the bkg (free parameters). Also add the rest of the BKG functions
1579     {
1580       fitType += Form(":kVWG=%f",minvResult->GetValue("kVWG"));
1581       fitType += Form(":mVWG=%f",minvResult->GetValue("mVWG"));
1582       fitType += Form(":sVWG1=%f",minvResult->GetValue("sVWG1"));
1583       fitType += Form(":sVWG2=%f",minvResult->GetValue("sVWG2"));
1584
1585       msg += " + VWG Bkg parameters";
1586     }
1587     else if ( fitType.Contains("POL2EXP_") || fitType.Contains("POL2EXPINDEPTAILS") )
1588     {
1589       fitType += Form(":kPol2Exp=%f",minvResult->GetValue("kPol2Exp"));
1590       fitType += Form(":pol0=%f",minvResult->GetValue("pol0"));
1591       fitType += Form(":pol1=%f",minvResult->GetValue("pol1"));
1592       fitType += Form(":pol2=%f",minvResult->GetValue("pol2"));
1593       fitType += Form(":exp=%f",minvResult->GetValue("exp"));
1594
1595       msg += " + Pol2xExp Bkg parameters";
1596     }
1597     else if ( fitType.Contains("POL4EXP_") || fitType.Contains("POL4EXPINDEPTAILS") )
1598     {
1599       fitType += Form(":kPol4Exp=%f",minvResult->GetValue("kPol4Exp"));
1600       fitType += Form(":pol0=%f",minvResult->GetValue("pol0"));
1601       fitType += Form(":pol1=%f",minvResult->GetValue("pol1"));
1602       fitType += Form(":pol2=%f",minvResult->GetValue("pol2"));
1603       fitType += Form(":pol3=%f",minvResult->GetValue("pol3"));
1604       fitType += Form(":pol4=%f",minvResult->GetValue("pol4"));
1605       fitType += Form(":exp=%f",minvResult->GetValue("exp"));
1606
1607       msg += " + Pol4xExp Bkg parameters";
1608     }
1609     std::cout << "Using " << msg.Data() << " from " << minvResult->GetName() <<  " inv mass result" << std::endl;
1610     std::cout << "" << std::endl;
1611   }
1612   else
1613   {
1614     AliError(Form("Cannot get tails from result. Result %s not found",minvResult->GetName()));
1615     fitType = "";
1616     return;
1617   }
1618 }
1619
1620 //_____________________________________________________________________________
1621 ULong64_t AliAnalysisMuMu::GetTriggerScalerCount(const char* triggerList, Int_t runNumber)
1622 {
1623   // Get the expected (from OCDB scalers) trigger count
1624
1625   AliAnalysisTriggerScalers ts(runNumber,Config()->OCDBPath());
1626
1627   TObjArray* triggers = TString(triggerList).Tokenize(",");
1628   TObjString* trigger;
1629   TIter next(triggers);
1630   ULong64_t n(0);
1631
1632   while ( ( trigger = static_cast<TObjString*>(next()) ) )
1633   {
1634     AliAnalysisTriggerScalerItem* item = ts.GetTriggerScaler(runNumber,"L2A",trigger->String().Data());
1635     if (item)
1636     {
1637       n += item->Value();
1638     }
1639     delete item;
1640   }
1641   delete triggers;
1642
1643   return n;
1644 }
1645
1646 //_____________________________________________________________________________
1647 AliAnalysisMuMuSpectra* AliAnalysisMuMu::GetSpectra(const char* what, const char* flavour) const
1648 {
1649   /// get a given spectra
1650
1651   TString swhat(what);
1652   TString sflavour(flavour);
1653   swhat.ToUpper();
1654   sflavour.ToUpper();
1655
1656   TString spectraName(Form("/%s/%s/PP/%s/PSI-%s",
1657                            First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,IsSimulation())).Data(),
1658                            First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,IsSimulation())).Data(),
1659                            First(Config()->GetList(AliAnalysisMuMuConfig::kPairSelectionList,IsSimulation())).Data(),
1660                            swhat.Data()));
1661
1662   if (sflavour.Length()>0)
1663   {
1664     spectraName += "-";
1665     spectraName += sflavour.Data();
1666   }
1667
1668   return SPECTRA(spectraName.Data());
1669 }
1670
1671 //_____________________________________________________________________________
1672 TH1* AliAnalysisMuMu::PlotAccEfficiency(const char* whatever)
1673 {
1674   //FIXME::Make it general
1675   if ( !IsSimulation() )
1676   {
1677     AliError("Could not get AccxEff histo: Not simulation file");
1678     return 0x0;
1679   }
1680
1681   TString path(Form("/%s/%s/%s/%s",
1682                     First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kTRUE)).Data(),
1683                     First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kTRUE)).Data(),
1684                     First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,kTRUE)).Data(),
1685                     First(Config()->GetList(AliAnalysisMuMuConfig::kPairSelectionList,kTRUE)).Data()));
1686
1687   AliAnalysisMuMuSpectra* s = static_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(Form("%s/%s",path.Data(),whatever)));
1688   if ( !s )
1689   {
1690     AliError(Form("No AccEff spectra %s found in %s",whatever,path.Data()));
1691     return 0x0;
1692   }
1693
1694   return s->Plot(Form("AccEff%s",GetParticleName()),"PSICOUNT",kFALSE);//_2.2_3.9
1695
1696 }
1697
1698 //_____________________________________________________________________________
1699 TH1* AliAnalysisMuMu::PlotSystematicsTestsRelative(const char* quantity,const char* flavour,const char* value2Test)
1700 {
1701   /// what,quantity and flavour defines de binning to test (output will be 1 plot per bin)
1702   /// value2test is the observable we want to test ( i.e. NJpsi(bin)/integratedNJpsi, <pt>(bin)/integrated<pt>... )
1703   /// --value2test == yield -> NJpsi(bin)/integratedNJpsi
1704   /// --value2test == mpt -> <pt>(bin)/integrated<pt>
1705
1706   TString path(Form("%s/%s/%s/%s",
1707                     First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kFALSE)).Data(),
1708                     First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data(),
1709                     First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,kFALSE)).Data(),
1710                     First(Config()->GetList(AliAnalysisMuMuConfig::kPairSelectionList,kFALSE)).Data()));
1711
1712   TString svalue2Test(value2Test);
1713   TString shName("");
1714   TString sObsInt("");
1715   TString sObsBin("");
1716   TString sflavour(flavour);
1717   TString squantity(quantity);
1718   squantity.ToUpper();
1719
1720   if ( !svalue2Test.CompareTo("yield",TString::kIgnoreCase) )
1721   {
1722     sObsInt = "PSI-INTEGRATED-AccEffCorr";
1723     sObsBin = Form("PSI-%s-AccEffCorr",squantity.Data());
1724     svalue2Test = "NofJPsi";
1725     shName = "N^{J/#psi}_{bin}/N^{J/#psi}_{int}";
1726   }
1727   else if ( !svalue2Test.CompareTo("mpt",TString::kIgnoreCase) )
1728   {
1729     sObsInt = "PSI-INTEGRATED-AccEffCorr-MeanPtVsMinvUS";
1730     sObsBin = Form("PSI-%s-AccEffCorr-MeanPtVsMinvUS",squantity.Data());
1731     svalue2Test = "MeanPtJPsi";
1732     shName = "<p_{t}>^{bin}/<p_{t}>^{int}";
1733   }
1734   else
1735   {
1736     AliError("unrecognized value to test");
1737     return  0x0;
1738   }
1739
1740   TString id(Form("/TESTSYST/%s",path.Data()));
1741
1742   //--Get the integrated results
1743   AliAnalysisMuMuSpectra* sInt = static_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(Form("/%s/%s",path.Data(),sObsInt.Data())));
1744   if ( !sInt )
1745   {
1746     AliError(Form("No integrated spectra %s found in %s",sObsInt.Data(),path.Data()));
1747     return 0x0;
1748   }
1749
1750   TObjArray* bin = BIN()->CreateBinObjArray("psi","integrated","");
1751   if ( !bin )
1752   {
1753     AliError("No integrated bin found");
1754     return 0x0;
1755   }
1756   AliAnalysisMuMuBinning::Range* r = static_cast<AliAnalysisMuMuBinning::Range*>(bin->At(0));
1757
1758   AliAnalysisMuMuJpsiResult* resInt =  static_cast<AliAnalysisMuMuJpsiResult*>(sInt->GetResultForBin(*r));
1759   if ( !resInt )
1760   {
1761     AliError(Form("No integrated result found in spectra %s at %s",sObsInt.Data(),path.Data()));
1762     return 0x0;
1763   }
1764   TObjArray* sresIntArray = resInt->SubResults();
1765 //  Int_t nsresInt = sresIntArray->GetEntriesFast();
1766
1767   bin = BIN()->CreateBinObjArray("psi",squantity.Data(),sflavour.Data());//memory leak?
1768   if ( !bin )
1769   {
1770     AliError(Form("%s-%s-%s binning does not exist","psi",squantity.Data(),sflavour.Data()));
1771     return 0x0;
1772   }
1773
1774   Int_t nbin = bin->GetEntries();
1775   Int_t nsres = sresIntArray->GetEntries();
1776   TObjArray* sResultNameArray= new TObjArray();
1777   std::vector<std::vector<double> > valuesArr;
1778   valuesArr.resize(nbin+1, std::vector<double>(nsres,0));
1779   std::vector<std::vector<double> > valuesErrorArr;
1780   valuesErrorArr.resize(nbin+1, std::vector<double>(nsres,0));
1781
1782   TIter nextIntSubResult(sresIntArray);
1783   AliAnalysisMuMuResult* sresInt(0x0);
1784   Int_t nBkgParametrizations(0);
1785   Int_t j(0);
1786   while ( ( sresInt = static_cast<AliAnalysisMuMuResult*>(nextIntSubResult()) ) ) //Integrated SubResult loop
1787   {
1788     TString sresIntName(sresInt->GetName());
1789
1790     valuesArr[0][j] = sresInt->GetValue(svalue2Test.Data());
1791     valuesErrorArr[0][j] = sresInt->GetErrorStat(svalue2Test.Data());
1792     sResultNameArray->Add(new TObjString(sresIntName.Data()));
1793
1794     if ( sresIntName.Contains("CB2") )
1795     {
1796       nBkgParametrizations++;
1797     }
1798     j++;
1799   }
1800
1801   //--Get the bin per bin results
1802   AliAnalysisMuMuSpectra* sBin = static_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(Form("/%s/%s",path.Data(),sObsBin.Data())));
1803   if ( !sBin )
1804   {
1805     AliError(Form("No integrated spectra %s found in %s",sObsBin.Data(),path.Data()));
1806     return 0x0;
1807     delete bin;
1808   }
1809
1810   TIter nextBin(bin);
1811   AliAnalysisMuMuJpsiResult* sresBin(0x0);
1812   Int_t i(1);
1813   while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) ) //Bin loop
1814   {
1815     sresBin =  static_cast<AliAnalysisMuMuJpsiResult*>(sBin->GetResultForBin(*r));
1816     if ( !sresBin )
1817     {
1818       AliError(Form("No result found in spectra %s at %s for bin %s",sObsInt.Data(),path.Data(),r->AsString().Data()));
1819       return 0x0;
1820       delete bin;
1821     }
1822
1823     TObjArray* sresBinArray = sresBin->SubResults();
1824
1825     TIter nextSubResult(sresBinArray);
1826     j = 0;
1827     while ( (sresBin = static_cast<AliAnalysisMuMuJpsiResult*>(nextSubResult())) ) // Subresults loop
1828     {
1829       valuesArr[i][j] = sresBin->GetValue(svalue2Test.Data());
1830       valuesErrorArr[i][j] = sresBin->GetErrorStat(svalue2Test.Data());
1831
1832       j++;
1833
1834     } //End Subresults loop
1835     i++;
1836
1837   } //End bin loop
1838
1839   TH1* hsyst = new TH1F(Form("%s_Systematics",value2Test),Form("%s Systematics results",value2Test),nbin,0,nbin);
1840
1841   TString binName("");
1842   for ( Int_t b = 1 ; b < i ; b++ ) //Bin loop
1843   {
1844     binName = static_cast<AliAnalysisMuMuBinning::Range*>(bin->At(b-1))->AsString().Data();
1845     TString savePath(Form("%s/%s",id.Data(),binName.Data()));
1846
1847     Int_t binUCTotal(1);
1848
1849     TH1* hratiosBin = new TH1D(Form("SystTests_%s_Bkg_%s",sObsBin.Data(),binName.Data()),
1850                           Form("%s Systematics tests for %s",binName.Data(),shName.Data()),j*nBkgParametrizations,0,
1851                           j*nBkgParametrizations);
1852
1853     for ( Int_t k = 0 ; k <= j ; k++ )
1854     {
1855       TString hName("");
1856
1857       if ( k == 0 ) hName = "UC";
1858       else hName = sResultNameArray->At(k-1)->GetName();
1859
1860
1861       TH1* hratios = new TH1D(Form("SystTests_%s_%s",sObsBin.Data(),hName.Data()),
1862                               Form("%s Systematics tests for %s(%s)",binName.Data(),shName.Data(),hName.Data()),j,0,j);
1863       hratios->SetNdivisions(j+1);
1864
1865       TH1* hratiosUC(0x0);
1866       if ( k != 0 )
1867       {
1868         hratiosUC = new TH1D(Form("SystTests_%s_%s_Bkg",sObsBin.Data(),hName.Data()),
1869                                          Form("%s Systematics tests for %s(%s bkg)",binName.Data(),shName.Data(),hName.Data()),nBkgParametrizations,0,nBkgParametrizations);
1870         hratiosUC->SetNdivisions(nBkgParametrizations+1);
1871       }
1872
1873       TString signalName(hName.Data());
1874       Int_t sizeName = signalName.Sizeof();
1875       signalName.Remove(2,sizeName-3);
1876       Int_t binUC(1);
1877
1878       for ( Int_t l = 0; l < j ; l++) //Subresults loop
1879       {
1880         TString binSignalName(sResultNameArray->At(l)->GetName());
1881
1882         Double_t ratio,ratioError;
1883         if ( k == 0)
1884         {
1885           ratio = valuesArr[b][l] / valuesArr[0][l];
1886           ratioError = TMath::Sqrt( TMath::Power(valuesErrorArr[b][l] / valuesArr[0][l],2.) + TMath::Power(valuesArr[b][l]*valuesErrorArr[0][l] / TMath::Power(valuesArr[0][l],2.),2.) );
1887
1888           hratios->GetXaxis()->SetBinLabel(l+1,sResultNameArray->At(l)->GetName());
1889         }
1890         else
1891         {
1892           ratio = valuesArr[b][l] / valuesArr[0][k-1];
1893           ratioError = TMath::Sqrt( TMath::Power(valuesErrorArr[b][l] / valuesArr[0][k-1],2.) + TMath::Power(valuesArr[b][l]*valuesErrorArr[0][k-1] / TMath::Power(valuesArr[0][k-1],2.),2.) );
1894
1895           hratios->GetXaxis()->SetBinLabel(l+1,Form("%s/%s",sResultNameArray->At(l)->GetName(),sResultNameArray->At(k-1)->GetName()));
1896
1897           if ( binSignalName.Contains(signalName.Data()) )
1898           {
1899             hratiosUC->GetXaxis()->SetBinLabel(binUC,Form("%s/%s",sResultNameArray->At(l)->GetName(),sResultNameArray->At(k-1)->GetName()));
1900
1901             hratiosUC->SetBinContent(binUC,ratio);
1902             hratiosUC->SetBinError(binUC,ratioError);
1903
1904             hratiosBin->GetXaxis()->SetBinLabel(binUCTotal,Form("%s/%s",sResultNameArray->At(l)->GetName(),sResultNameArray->At(k-1)->GetName()));
1905
1906             hratiosBin->SetBinContent(binUCTotal,ratio);
1907             hratiosBin->SetBinError(binUCTotal,ratioError);
1908
1909             binUC++;
1910             binUCTotal++;
1911           }
1912         }
1913
1914         hratios->SetBinContent(l+1,ratio);
1915         hratios->SetBinError(l+1,ratioError);
1916       }
1917
1918
1919
1920
1921 //      TH1* o = OC()->Histo(Form("%s",savePath.Data()),hratios->GetName());
1922 //
1923 //      if (o)
1924 //      {
1925 //        AliWarning(Form("Replacing %s/%s",savePath.Data(),hratios->GetName()));
1926 //        OC()->Remove(Form("%s/%s",savePath.Data(),hratios->GetName()));
1927 //      }
1928 //
1929 //      Bool_t adoptOK = OC()->Adopt(savePath.Data(),hratios);
1930 //
1931 //      if ( adoptOK ) std::cout << "+++syst histo " << hratios->GetName() << " adopted" << std::endl;
1932 //      else AliError(Form("Could not adopt syst histo %s",hratios->GetName()));
1933 //
1934 //      if ( hratiosUC )
1935 //      {
1936 //        o = OC()->Histo(Form("%s",savePath.Data()),hratiosUC->GetName());
1937 //
1938 //        if (o)
1939 //        {
1940 //          AliWarning(Form("Replacing %s/%s",savePath.Data(),hratiosUC->GetName()));
1941 //          OC()->Remove(Form("%s/%s",savePath.Data(),hratiosUC->GetName()));
1942 //        }
1943 //
1944 //        adoptOK = OC()->Adopt(savePath.Data(),hratiosUC);
1945 //
1946 //        if ( adoptOK ) std::cout << "+++syst histo " << hratiosUC->GetName() << " adopted" << std::endl;
1947 //        else AliError(Form("Could not adopt syst histo %s",hratiosUC->GetName()));
1948 //      }
1949     }
1950
1951
1952     //Syst computation for bin
1953     Double_t num(0.),deno(0.);
1954     for ( Int_t m = 1 ; m <= hratiosBin->GetNbinsX() ; m++ )
1955     {
1956       Double_t value = hratiosBin->GetBinContent(m);
1957       Double_t error = hratiosBin->GetBinError(m);
1958
1959       num += value/TMath::Power(error,2.);
1960       deno += 1./TMath::Power(error,2.);
1961     }
1962
1963     Double_t mean = num/deno;
1964     Double_t v1(0.),v2(0.),sum(0.);
1965     for ( Int_t l = 1 ; l <= hratiosBin->GetNbinsX() ; l++ )
1966     {
1967       Double_t value = hratiosBin->GetBinContent(l);
1968       Double_t error = hratiosBin->GetBinError(l);
1969
1970       Double_t wi = 1./TMath::Power(error,2.);
1971       v1 += wi;
1972       v2 += wi*wi;
1973       Double_t diff = value - mean;
1974       sum += wi*diff*diff;
1975
1976     }
1977
1978     Double_t syst = TMath::Sqrt( (v1/(v1*v1-v2)) * sum);
1979
1980     hsyst->GetXaxis()->SetBinLabel(b,binName.Data());
1981     hsyst->SetBinContent(b,(syst*100.)/mean);
1982
1983
1984     //___
1985     TF1* meanF = new TF1("mean","[0]",0,j*nBkgParametrizations);
1986     meanF->SetParameter(0,mean);
1987
1988     TF1* meanFPS = new TF1("meanPS","[0]",0,j*nBkgParametrizations);
1989     meanFPS->SetParameter(0,mean+syst);
1990     meanFPS->SetLineStyle(2);
1991
1992     TF1* meanFMS = new TF1("meanMS","[0]",0,j*nBkgParametrizations);
1993     meanFMS->SetParameter(0,mean-syst);
1994     meanFMS->SetLineStyle(2);
1995
1996     hratiosBin->GetListOfFunctions()->Add(meanF);
1997     hratiosBin->GetListOfFunctions()->Add(meanFPS);
1998     hratiosBin->GetListOfFunctions()->Add(meanFMS);
1999
2000     TH1* o = OC()->Histo(Form("%s",savePath.Data()),hratiosBin->GetName());
2001
2002     if (o)
2003     {
2004       AliWarning(Form("Replacing %s/%s",savePath.Data(),hratiosBin->GetName()));
2005       OC()->Remove(Form("%s/%s",savePath.Data(),hratiosBin->GetName()));
2006     }
2007
2008     Bool_t adoptOK = OC()->Adopt(savePath.Data(),hratiosBin);
2009
2010     if ( adoptOK ) std::cout << "+++syst histo " << hratiosBin->GetName() << " adopted" << std::endl;
2011     else AliError(Form("Could not adopt syst histo %s",hratiosBin->GetName()));
2012   }
2013
2014
2015   TH1* o = OC()->Histo(Form("%s",id.Data()),hsyst->GetName());
2016
2017   if (o)
2018   {
2019     AliWarning(Form("Replacing %s/%s",id.Data(),hsyst->GetName()));
2020     OC()->Remove(Form("%s/%s",id.Data(),hsyst->GetName()));
2021   }
2022
2023   Bool_t adoptOK = OC()->Adopt(id.Data(),hsyst);
2024
2025   if ( adoptOK ) std::cout << "+++syst histo " << hsyst->GetName() << " adopted" << std::endl;
2026   else AliError(Form("Could not adopt syst histo %s",hsyst->GetName()));
2027
2028
2029   delete bin;
2030   delete sResultNameArray;
2031
2032   return 0x0;
2033
2034 }
2035
2036
2037 //_____________________________________________________________________________
2038 TH1* AliAnalysisMuMu::PlotJpsiYield(const char* whatever)
2039 {
2040
2041   //FIXME::Make it general
2042   if ( IsSimulation() )
2043   {
2044     AliError("Cannot compute J/Psi yield: Is a simulation file");
2045     return 0x0;
2046   }
2047
2048   TString path(Form("/%s/%s/%s/%s",
2049                     First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kFALSE)).Data(),
2050                     First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data(),
2051                     First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,kFALSE)).Data(),
2052                     First(Config()->GetList(AliAnalysisMuMuConfig::kPairSelectionList,kFALSE)).Data()));
2053
2054   AliAnalysisMuMuSpectra* s = static_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(Form("%s/%s",path.Data(),whatever)));
2055   if ( !s )
2056   {
2057     AliError(Form("No spectra %s found in %s",whatever,path.Data()));
2058     return 0x0;
2059   }
2060
2061   std::cout << "Number of J/Psi:" << std::endl;
2062   TH1* hry = s->Plot("NofJPsi","PSIPSIPRIMECB2VWGINDEPTAILS",kFALSE); //Number of Jpsi
2063   std::cout << "" << std::endl;
2064
2065   std::cout << "Equivalent number of MB events:" << std::endl;
2066   TH1* hN = ComputeEquNofMB();
2067   std::cout << "" << std::endl;
2068
2069   TH1* hy = static_cast<TH1*>(hry->Clone("CorrJPsiYields"));
2070   Double_t bR = 0.0593; // BR(JPsi->mu+mu-)
2071   Double_t bRerror = 0.0006 ;
2072
2073   for (Int_t i = 1 ; i <= hy->GetNbinsX() ; i++)
2074   {
2075     Double_t yield = hry->GetBinContent(i)/(hN->GetBinContent(i)*bR);
2076     Double_t yieldError = TMath::Sqrt(TMath::Power(hry->GetBinError(i)/(hN->GetBinContent(i)*bR),2.) +
2077                                       TMath::Power(hN->GetBinError(i)*bR/TMath::Power(hN->GetBinContent(i)*bR,2.),2.) +
2078                                       TMath::Power(hry->GetBinContent(i)*hN->GetBinContent(i)*bRerror/TMath::Power(hN->GetBinContent(i)*bR,2.),2.));
2079
2080     std::cout << yield << " +- " << yieldError << std::endl;
2081
2082     hy->SetBinContent(i,yield);
2083     hy->SetBinError(i,yieldError);
2084   }
2085
2086   delete hry;
2087   delete hN;
2088
2089   return hy;
2090 }
2091
2092
2093 //_____________________________________________________________________________
2094 UInt_t AliAnalysisMuMu::GetSum(AliCounterCollection& cc, const char* triggerList,
2095                                const char* eventSelection, Int_t runNumber)
2096 {
2097   TObjArray* ktrigger = cc.GetKeyWords("trigger").Tokenize(",");
2098   TObjArray* kevent = cc.GetKeyWords("event").Tokenize(",");
2099   TObjArray* a = TString(triggerList).Tokenize(" ");
2100   TIter next(a);
2101   TObjString* str;
2102
2103   UInt_t n(0);
2104
2105   TString sEventSelection(eventSelection);
2106   sEventSelection.ToUpper();
2107
2108   if ( kevent->FindObject(sEventSelection.Data()) )
2109   {
2110     while ( ( str = static_cast<TObjString*>(next()) ) )
2111     {
2112       if ( ktrigger->FindObject(str->String().Data()) )
2113       {
2114         if ( runNumber < 0 )
2115         {
2116           n +=  static_cast<UInt_t>(cc.GetSum(Form("trigger:%s/event:%s",str->String().Data(),eventSelection)));
2117         }
2118         else
2119         {
2120           n +=  static_cast<UInt_t>(cc.GetSum(Form("trigger:%s/event:%s/run:%d",str->String().Data(),eventSelection,runNumber)));
2121         }
2122       }
2123     }
2124   }
2125
2126   delete a;
2127   delete ktrigger;
2128   delete kevent;
2129   return n;
2130 }
2131
2132 //_____________________________________________________________________________
2133 Bool_t
2134 AliAnalysisMuMu::GetCollections(const char* rootfile,
2135                                 AliMergeableCollection*& oc,
2136                                 AliCounterCollection*& cc,
2137                                 AliAnalysisMuMuBinning*& bin,
2138                                 std::set<int>& runnumbers)
2139 {
2140   oc = 0x0;
2141   cc = 0x0;
2142   bin = 0x0;
2143
2144   TFile* f = static_cast<TFile*>(gROOT->GetListOfFiles()->FindObject(rootfile));
2145   if (!f)
2146   {
2147     f = TFile::Open(rootfile);
2148   }
2149
2150   if ( !f || f->IsZombie() )
2151   {
2152     return kFALSE;
2153   }
2154
2155   f->GetObject("OC",oc);
2156   if (!oc)
2157   {
2158     f->GetObject("MC",oc);
2159   }
2160   f->GetObject("CC",cc);
2161
2162   TIter next(f->GetListOfKeys());
2163   TKey* key;
2164
2165   while ( ( key = static_cast<TKey*>(next())) && !bin )
2166   {
2167     if ( strcmp(key->GetClassName(),"AliAnalysisMuMuBinning")==0 )
2168     {
2169       bin = dynamic_cast<AliAnalysisMuMuBinning*>(key->ReadObj());
2170     }
2171   }
2172
2173   if (!oc || !cc)
2174   {
2175     AliErrorClass("Old file. Please upgrade it!");
2176
2177     return kFALSE;
2178   }
2179
2180   // get run list
2181   TObjArray* runs = cc->GetKeyWords("run").Tokenize(",");
2182   runs->Sort();
2183   TIter nextRun(runs);
2184   TObjString* srun;
2185
2186   runnumbers.clear();
2187
2188   while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
2189   {
2190     runnumbers.insert(srun->String().Atoi());
2191   }
2192
2193   delete runs;
2194
2195   return kTRUE;
2196 }
2197
2198 //_____________________________________________________________________________
2199 Bool_t AliAnalysisMuMu::IsSimulation() const
2200 {
2201   // whether or not we have MC information
2202
2203 //  return kFALSE;
2204
2205   if (!fMergeableCollection) return kFALSE;
2206
2207   TList* list = fMergeableCollection->CreateListOfKeys(0);
2208   TIter next(list);
2209   TObjString* str;
2210   Bool_t ok(kFALSE);
2211
2212   while ( ( str = static_cast<TObjString*>(next()) ) )
2213   {
2214     if ( str->String().Contains(AliAnalysisMuMuBase::MCInputPrefix()) ) ok = kTRUE;
2215   }
2216   delete list;
2217
2218   return ok;
2219 }
2220
2221 //_____________________________________________________________________________
2222 Int_t
2223 AliAnalysisMuMu::Jpsi(const char* what, const char* binningFlavour, Bool_t fitmPt)
2224 {
2225   // Fit the J/psi (and psiprime) peaks for the triggers in fDimuonTriggers list
2226   // what="integrated" => fit only fully integrated MinvUS
2227   // what="pt" => fit MinvUS in pt bins
2228   // what="y" => fit MinvUS in y bins
2229   // what="pt,y" => fit MinvUS in (pt,y) bins
2230
2231   TStopwatch timer;
2232
2233   if (!fMergeableCollection)
2234   {
2235     AliError("No mergeable collection. Consider Upgrade()");
2236     return 0;
2237   }
2238
2239   Int_t nfits(0);
2240
2241   TObjArray* triggerArray = Config()->GetListElements(AliAnalysisMuMuConfig::kDimuonTriggerList,IsSimulation());
2242   TObjArray* eventTypeArray = Config()->GetListElements(AliAnalysisMuMuConfig::kEventSelectionList,IsSimulation());
2243   TObjArray* pairCutArray = Config()->GetListElements(AliAnalysisMuMuConfig::kPairSelectionList,IsSimulation());
2244   TObjArray* whatArray = TString(what).Tokenize(",");
2245   TObjArray* centralityArray = Config()->GetListElements(AliAnalysisMuMuConfig::kCentralitySelectionList,IsSimulation());
2246
2247   TIter nextTrigger(triggerArray);
2248   TIter nextEventType(eventTypeArray);
2249   TIter nextPairCut(pairCutArray);
2250   TIter nextWhat(whatArray);
2251   TIter nextCentrality(centralityArray);
2252
2253   TObjString* trigger;
2254   TObjString* eventType;
2255   TObjString* pairCut;
2256   TObjString* swhat;
2257   TObjString* centrality;
2258
2259   while ( ( swhat = static_cast<TObjString*>(nextWhat()) ) )
2260   {
2261     AliAnalysisMuMuBinning* binning(0x0);
2262
2263     if ( fBinning && swhat->String().Length() > 0 )
2264     {
2265       binning = fBinning->Project("psi",swhat->String().Data(),binningFlavour);
2266     }
2267     else
2268     {
2269       binning = new AliAnalysisMuMuBinning;
2270       binning->AddBin("psi",swhat->String().Data());
2271     }
2272
2273     StdoutToAliDebug(1,std::cout << "++++++++++++ swhat=" << swhat->String().Data() << std::endl;);
2274
2275     std::cout << "" << std::endl;
2276     std::cout << "++++++++++++++++++" << "NEW BIN TYPE" << "+++++++++++++++++++" << std::endl;
2277     std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
2278     std::cout << "++++++++++++ swhat=" << swhat->String().Data() << "++++++++++++++++++++" << std::endl;
2279     std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
2280     std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
2281     std::cout << "" << std::endl;
2282     std::cout << "" << std::endl;
2283
2284     if (!binning)
2285     {
2286       AliError("oups. binning is NULL");
2287       continue;
2288     }
2289
2290     StdoutToAliDebug(1,binning->Print(););
2291
2292     nextTrigger.Reset();
2293
2294     while ( ( trigger = static_cast<TObjString*>(nextTrigger())) )
2295     {
2296       AliDebug(1,Form("TRIGGER %s",trigger->String().Data()));
2297
2298       nextEventType.Reset();
2299
2300       while ( ( eventType = static_cast<TObjString*>(nextEventType())) )
2301       {
2302         AliDebug(1,Form("--EVENTTYPE %s",eventType->String().Data()));
2303
2304         nextPairCut.Reset();
2305
2306         while ( ( pairCut = static_cast<TObjString*>(nextPairCut())) )
2307         {
2308           AliDebug(1,Form("----PAIRCUT %s",pairCut->String().Data()));
2309
2310           nextCentrality.Reset();
2311
2312           while ( ( centrality = static_cast<TObjString*>(nextCentrality()) ) )
2313           {
2314             AliDebug(1,"----Fitting...");
2315
2316             AliAnalysisMuMuSpectra* spectra = FitParticle("psi",
2317                                                           trigger->String().Data(), //Uncomment
2318                                                           eventType->String().Data(),
2319                                                           pairCut->String().Data(),
2320                                                           centrality->String().Data(),
2321                                                           *binning);
2322
2323             AliDebug(1,Form("----fitting done spectra = %p",spectra));
2324
2325             TString id(Form("/%s/%s/%s/%s",eventType->String().Data(),
2326                             trigger->String().Data(),
2327                             centrality->String().Data(),
2328                             pairCut->String().Data()));
2329             TObject* o;
2330
2331             if ( spectra )
2332             {
2333               ++nfits;
2334
2335               o = fMergeableCollection->GetObject(id.Data(),spectra->GetName());
2336
2337               AliDebug(1,Form("----nfits=%d id=%s o=%p",nfits,id.Data(),o));
2338
2339               if (o)
2340               {
2341                 AliWarning(Form("Replacing %s/%s",id.Data(),spectra->GetName()));
2342                 fMergeableCollection->Remove(Form("%s/%s",id.Data(),spectra->GetName()));
2343               }
2344
2345               Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),spectra);
2346
2347               if ( adoptOK ) std::cout << "+++Spectra " << spectra->GetName() << " adopted" << std::endl;
2348               else AliError(Form("Could not adopt spectra %s",spectra->GetName()));
2349
2350               StdoutToAliDebug(1,spectra->Print(););
2351             }
2352             else AliError("Error creating spectra"); //Uncomment
2353
2354             AliDebug(1,"----Fitting corrected spectra...");
2355
2356             AliAnalysisMuMuSpectra* spectraCorr = FitParticle("psi",
2357                                                               trigger->String().Data(),
2358                                                               eventType->String().Data(),
2359                                                               pairCut->String().Data(),
2360                                                               centrality->String().Data(),
2361                                                               *binning,"minv",kTRUE);
2362
2363             AliDebug(1,Form("----fitting done corrected spectra = %p",spectraCorr));
2364
2365             o = 0x0;
2366             if ( spectraCorr )
2367             {
2368               ++nfits;
2369
2370               o = fMergeableCollection->GetObject(id.Data(),spectraCorr->GetName());
2371
2372               AliDebug(1,Form("----nfits=%d id=%s o=%p",nfits,id.Data(),o));
2373
2374               if (o)
2375               {
2376                 AliWarning(Form("Replacing %s/%s",id.Data(),spectraCorr->GetName()));
2377                 fMergeableCollection->Remove(Form("%s/%s",id.Data(),spectraCorr->GetName()));
2378               }
2379
2380               Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),spectraCorr);
2381
2382               if ( adoptOK ) std::cout << "+++Spectra " << spectraCorr->GetName() << " adopted" << std::endl;
2383               else AliError(Form("Could not adopt spectra %s",spectraCorr->GetName()));
2384
2385               StdoutToAliDebug(1,spectraCorr->Print(););
2386             }
2387             else AliError("Error creating spectra");
2388
2389
2390             if (fitmPt)
2391             {
2392               AliDebug(1,"----Fitting mean pt...");
2393
2394               std::cout << "" << std::endl;
2395               std::cout << "" << std::endl;
2396               std::cout << "++++++++++++ Fitting mean Pt for " << swhat->String().Data() << " " << "slices" << std::endl; //Uncomment
2397
2398               if ( spectra )
2399               {
2400                 AliAnalysisMuMuSpectra* spectraMeanPt = FitParticle("psi",
2401                                                                   trigger->String().Data(),
2402                                                                   eventType->String().Data(),
2403                                                                   pairCut->String().Data(),
2404                                                                   centrality->String().Data(),
2405                                                                   *binning,"mpt"/*,*spectra*/);
2406
2407
2408
2409                 AliDebug(1,Form("----fitting done spectra = %p",spectraMeanPt));
2410                 o = 0x0;
2411
2412                 if ( spectraMeanPt )
2413                 {
2414                   ++nfits; //Review this
2415
2416                   o = fMergeableCollection->GetObject(id.Data(),spectraMeanPt->GetName());
2417
2418                   AliDebug(1,Form("----nfits=%d id=%s o=%p",nfits,id.Data(),o));
2419
2420                   if (o)
2421                   {
2422                     AliWarning(Form("Replacing %s/%s",id.Data(),spectraMeanPt->GetName()));
2423                     fMergeableCollection->Remove(Form("%s/%s",id.Data(),spectraMeanPt->GetName()));
2424                   }
2425
2426                   Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),spectraMeanPt);
2427                   //                  spectraMeanPt->Print();
2428                   if ( adoptOK ) std::cout << "+++Spectra " << spectraMeanPt->GetName() << " adopted" << std::endl;
2429                   else AliError(Form("Could not adopt spectra %s",spectraMeanPt->GetName()));
2430                 }
2431                 else AliError("Error creating spectra");
2432
2433               }
2434               else std::cout << "Mean pt fit failed: No inv mass spectra for " << swhat->String().Data() << " " << "slices" << std::endl; //Uncomment
2435
2436
2437               std::cout << "++++++++++++ Fitting corrected mean Pt for" << " " << swhat->String().Data() << " " << "slices" << std::endl;
2438
2439               if ( spectraCorr )
2440               {
2441                 AliAnalysisMuMuSpectra* spectraMeanPtCorr  = FitParticle("psi",
2442                                                                        trigger->String().Data(),
2443                                                                        eventType->String().Data(),
2444                                                                        pairCut->String().Data(),
2445                                                                        centrality->String().Data(),
2446                                                                        *binning,"mpt"/*,*spectraCorr*/,kTRUE);
2447
2448
2449
2450                 AliDebug(1,Form("----fitting done spectra = %p",spectraMeanPtCorr));
2451
2452                 o = 0x0;
2453
2454                 if ( spectraMeanPtCorr )
2455                 {
2456                   ++nfits; //Review this
2457
2458                   o = fMergeableCollection->GetObject(id.Data(),spectraMeanPtCorr->GetName());
2459
2460                   AliDebug(1,Form("----nfits=%d id=%s o=%p",nfits,id.Data(),o));
2461
2462                   if (o)
2463                   {
2464                     AliWarning(Form("Replacing %s/%s",id.Data(),spectraMeanPtCorr->GetName()));
2465                     fMergeableCollection->Remove(Form("%s/%s",id.Data(),spectraMeanPtCorr->GetName()));
2466                   }
2467
2468                   Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),spectraMeanPtCorr);
2469                   //                  spectraMeanPtCorr->Print();
2470                   if ( adoptOK ) std::cout << "+++Spectra " << spectraMeanPtCorr->GetName() << " adopted" << std::endl;
2471                   else AliError(Form("Could not adopt spectra %s",spectraMeanPtCorr->GetName()));
2472
2473                 }
2474                 else AliError("Error creating spectra");
2475
2476               }
2477
2478               else std::cout << "Corrected mean pt fit failed: No corrected inv mass spectra for " << swhat->String().Data() << " " << "slices" << std::endl;
2479
2480             }
2481           }
2482         }
2483       }
2484     }
2485   }
2486
2487   delete whatArray;
2488   delete triggerArray;
2489   delete eventTypeArray;
2490   delete pairCutArray;
2491   delete centralityArray;
2492
2493   StdoutToAliDebug(1,timer.Print(););
2494
2495   if (nfits)
2496   {
2497     Update();
2498 //    ReOpen(fFilename,"UPDATE");
2499 //    fMergeableCollection->Write("MC",TObjArray::kOverwrite);// | TObjArray::kSingleKey);
2500 //    ReOpen(fFilename,"READ");
2501   }
2502
2503
2504   return nfits;
2505
2506 }
2507
2508 //_____________________________________________________________________________
2509 TGraph* AliAnalysisMuMu::PlotEventSelectionEvolution(const char* trigger1, const char* event1,
2510                                                      const char* trigger2, const char* event2,
2511                                                      Bool_t drawFills,
2512                                                      Bool_t asRejection) const
2513 {
2514   if (!CC()) return 0x0;
2515
2516   const std::set<int>& runnumbers = RunNumbers();
2517
2518   TGraphErrors* g = new TGraphErrors(runnumbers.size());
2519
2520   std::set<int>::const_iterator it;
2521   Int_t i(0);
2522
2523   Double_t ymin(TMath::Limits<double>::Max());
2524   Double_t ymax(TMath::Limits<double>::Min());
2525
2526   for ( it = runnumbers.begin(); it != runnumbers.end(); ++it )
2527   {
2528     Int_t runNumber = *it;
2529     Double_t n = CC()->GetSum(Form("trigger:%s/event:%s/run:%d",trigger1,event1,runNumber));
2530     Double_t d = CC()->GetSum(Form("trigger:%s/event:%s/run:%d",trigger2,event2,runNumber));
2531     if (n>0 && d>0)
2532     {
2533       Double_t y = n/d;
2534
2535       if ( fCorrectionPerRun )
2536       {
2537         Double_t xcorr,ycorr;
2538         fCorrectionPerRun->GetPoint(i,xcorr,ycorr); // note that the fact that xcorr==runNumber has been checked by the SetCorrectionPerRun method
2539         y *= ycorr;
2540         // FIXME: should get the correction error here
2541       }
2542
2543       if ( asRejection ) y = 100*(1.0 - y);
2544       ymin = TMath::Min(ymin,y);
2545       ymax = TMath::Max(ymax,y);
2546       Double_t yerr = y*AliAnalysisMuMuResult::ErrorAB(n,TMath::Sqrt(n),d,TMath::Sqrt(d));
2547       g->SetPoint(i,runNumber,y);
2548       g->SetPointError(i,0.5,yerr);
2549
2550       ++i;
2551     }
2552
2553   }
2554
2555   TH2* hframe = new TH2F(Form("%s %s-%s / %s-%s",(asRejection ? "1 - ":""),trigger1,event1,trigger2,event2),
2556                          Form("%s %s-%s / %s-%s",(asRejection ? "1 - ":""),trigger1,event1,trigger2,event2),
2557                          runnumbers.size()+50,
2558                          *(runnumbers.begin())-25,
2559                          *(runnumbers.rbegin())+25,100,0,ymax*1.3);
2560
2561   gStyle->SetOptStat(0);
2562
2563   hframe->Draw();
2564
2565   hframe->GetXaxis()->SetNoExponent();
2566
2567   hframe->GetYaxis()->SetTitle(asRejection ? "Rejection (%)" : "Ratio");
2568
2569   g->Set(i);
2570   g->SetTitle(Form("%s %s-%s / %s-%s",(asRejection ? "1 - ":""),trigger1,event1,trigger2,event2));
2571   g->GetXaxis()->SetNoExponent();
2572   g->Draw("lp");
2573
2574   AliAnalysisTriggerScalers ts(RunNumbers(),Config()->OCDBPath());
2575
2576   if ( drawFills )
2577   {
2578     ts.DrawFills(ymin,ymax);
2579     g->Draw("lp");
2580   }
2581
2582
2583   std::map<std::string, std::pair<int,int> > periods;
2584
2585   ts.GetLHCPeriodBoundaries(periods);
2586
2587   TLegend* legend = new TLegend(0.15,0.82,0.90,0.92);
2588   legend->SetFillColor(0);
2589   Int_t n(0);
2590
2591
2592   for ( std::map<std::string, std::pair<int,int> >::const_iterator pit = periods.begin(); pit != periods.end(); ++pit )
2593   {
2594     std::string period = pit->first;
2595     int run1 = (pit->second).first;
2596     int run2 = (pit->second).second;
2597     int nruns(0);
2598     for ( std::set<int>::const_iterator rit = RunNumbers().begin(); rit != RunNumbers().end(); ++ rit )
2599     {
2600       if ( (*rit) >= run1 && (*rit) <= run2 )
2601       {
2602         ++nruns;
2603       }
2604     }
2605     AliInfo(Form("Period %s runs %6d-%6d ; %d actual runs",period.c_str(),run1,run2,nruns));
2606
2607     g->Fit("pol0","+Q","",run1,run2);
2608     TF1* func = static_cast<TF1*>(g->GetListOfFunctions()->Last());
2609     if (func)
2610     {
2611       func->SetLineColor(2+n);
2612       legend->AddEntry(func,Form("%s %5.2f #pm %5.2f %s (rel. error %5.2f %%)",period.c_str(),func->GetParameter(0),func->GetParError(0),
2613                                  (asRejection ? "%":""),100*func->GetParError(0)/func->GetParameter(0)));
2614       ++n;
2615     }
2616   }
2617
2618   legend->SetNColumns(3);
2619
2620   Double_t mean = TMath::Mean(g->GetN(),g->GetY());
2621   Double_t rms = TMath::RMS(g->GetN(),g->GetY());
2622
2623   legend->AddEntry("",Form("Mean %5.2f RMS %5.2f (%5.2f %%)",mean,rms,(mean) ? 100.0*rms/mean : 0.0),"");
2624
2625   legend->Draw();
2626
2627   return g;
2628 }
2629
2630 //_____________________________________________________________________________
2631 void AliAnalysisMuMu::ShowList(const char* title, const TString& list, const char separator) const
2632 {
2633   /// Show a list of strings
2634   TObjArray* parts = list.Tokenize(separator);
2635
2636   std::cout << title << " (" << parts->GetEntries() << ") : " << std::endl;
2637
2638   TIter next(parts);
2639   TObjString* str;
2640
2641   while ( ( str = static_cast<TObjString*>(next()) ) )
2642   {
2643     std::cout << "    " << str->String().Data() << std::endl;
2644   }
2645
2646   if ( parts->GetEntries()==0) std::cout << endl;
2647
2648   delete parts;
2649 }
2650
2651 //_____________________________________________________________________________
2652 void AliAnalysisMuMu::Print(Option_t* opt) const
2653 {
2654     /// printout
2655   std::cout << "Reading from file : " << fFilename.Data() << std::endl;
2656
2657   TString copt(opt);
2658
2659   if (IsSimulation() || SIM() )
2660   {
2661     copt += "SIM";
2662   }
2663
2664   if ( !IsSimulation() )
2665   {
2666     copt += "REAL";
2667   }
2668
2669   Config()->Print(copt.Data());
2670
2671   if ( RunNumbers().size() > 1 )
2672   {
2673     std::cout << RunNumbers().size() << " runs";
2674   }
2675   else
2676   {
2677     std::cout << RunNumbers().size() << " run";
2678   }
2679
2680   if ( fCorrectionPerRun )
2681   {
2682     std::cout << " with correction factors";
2683   }
2684   std::cout << std::endl;
2685   Int_t i(0);
2686   for ( std::set<int>::const_iterator it = RunNumbers().begin(); it != RunNumbers().end(); ++it )
2687   {
2688     std::cout << (*it);
2689     if ( fCorrectionPerRun )
2690     {
2691       std::cout << Form("(%e)",fCorrectionPerRun->GetY()[i]);
2692     }
2693     std::cout << ",";
2694     ++i;
2695   }
2696   std::cout << std::endl;
2697
2698   TString sopt(opt);
2699   sopt.ToUpper();
2700
2701   if ( sopt.Contains("BIN") && BIN() )
2702   {
2703     std::cout << "Binning : " << std::endl;
2704     TString topt(sopt);
2705     topt.ReplaceAll("BIN","");
2706     BIN()->Print(topt.Data());
2707   }
2708   if ( sopt.Contains("MC") && OC() )
2709   {
2710     TString topt(sopt);
2711     topt.ReplaceAll("MC","");
2712     OC()->Print(topt.Data());
2713   }
2714   if ( sopt.Contains("CC") && CC() )
2715   {
2716     CC()->Print("trigger/event");
2717   }
2718
2719   if ( sopt.Contains("SIZE") )
2720   {
2721     TFile* f = ReOpen(fFilename,"READ");
2722     TIter next(f->GetListOfKeys());
2723     TKey* key;
2724
2725     while ( ( key = static_cast<TKey*>(next()) ) )
2726     {
2727       std::cout << key->GetName() << " " << key->GetNbytes() << " " << key->GetObjlen() << std::endl;
2728     }
2729   }
2730 }
2731
2732 //_____________________________________________________________________________
2733 TFile* AliAnalysisMuMu::ReOpen(const char* filename, const char* mode) const
2734 {
2735   /// Tries to reopen the file with a new mode
2736
2737   TFile* f = static_cast<TFile*>(gROOT->GetListOfFiles()->FindObject(filename));
2738
2739   if (f)
2740   {
2741     delete f;
2742   }
2743
2744   f = TFile::Open(filename,mode);
2745
2746   if ( !f || !f->IsOpen() )
2747   {
2748     AliError(Form("Cannot open file %s in mode %s",filename,mode));
2749     return 0x0;
2750   }
2751
2752   return f;
2753 }
2754
2755 //_____________________________________________________________________________
2756 void AliAnalysisMuMu::SetCentralitySelectionList(const char* centralitySelectionList)
2757 {
2758   /// Set centralities to be used during fitting
2759   /// centralitySelectionList is a regular expression.
2760
2761   TObjArray* centralities = BIN()->CreateBinObjArray("centrality");
2762   TIter next(centralities);
2763   AliAnalysisMuMuBinning::Range* r;
2764
2765   TString csl;
2766
2767   TPRegexp re(centralitySelectionList);
2768
2769   while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(next()) ) )
2770   {
2771     AliDebug(1,Form("r=%s",r->AsString().Data()));
2772     if ( re.MatchB(r->AsString()) )
2773     {
2774       csl += r->AsString();
2775       csl += ",";
2776     }
2777   }
2778
2779   Config()->SetList(AliAnalysisMuMuConfig::kCentralitySelectionList,IsSimulation(),csl);
2780
2781   delete centralities;
2782 }
2783
2784 //_____________________________________________________________________________
2785 Bool_t AliAnalysisMuMu::SetCorrectionPerRun(const TGraph& corr, const char* formula)
2786 {
2787     /// Sets the graph used to correct values per run
2788   delete fCorrectionPerRun;
2789   fCorrectionPerRun=0x0;
2790
2791   // check that corr has the same runs as we do
2792
2793   Int_t i(0);
2794
2795   for ( std::set<int>::const_iterator it = RunNumbers().begin(); it != RunNumbers().end(); ++it )
2796   {
2797     Int_t corrRun = TMath::Nint(corr.GetX()[i]);
2798     if (corrRun != *it)
2799     {
2800       AliError(Form("%d-th run mistmatch %d vs %d",i,corrRun,*it));
2801
2802       return kFALSE;
2803     }
2804     ++i;
2805   }
2806
2807   fCorrectionPerRun = new TGraphErrors(corr.GetN());
2808
2809   TFormula* tformula(0x0);
2810   if ( strlen(formula) > 0 )
2811   {
2812     tformula = new TFormula("SetCorrectionPerRunFormula",formula);
2813   }
2814
2815   i = 0;
2816
2817   for ( std::set<int>::const_iterator it = RunNumbers().begin(); it != RunNumbers().end(); ++it )
2818   {
2819     Double_t y = corr.GetY()[i];
2820
2821     if ( tformula )
2822     {
2823       y = tformula->Eval(y);
2824     }
2825     fCorrectionPerRun->SetPoint(i,corr.GetX()[i],y);
2826     ++i;
2827   }
2828
2829   delete formula;
2830
2831   return kTRUE;
2832 }
2833
2834 //_____________________________________________________________________________
2835 void AliAnalysisMuMu::SetNofInputParticles(AliAnalysisMuMuJpsiResult& r)
2836 {
2837   /// Set the "NofInput" variable(s) of one result
2838
2839   TString hname(Form("MinvUS%s",r.Bin().AsString().Data()));
2840
2841   TH1* hinput = fMergeableCollection->Histo(Form("/%s/ALL/ANY/V0A/INYRANGE",AliAnalysisMuMuBase::MCInputPrefix()),hname.Data());
2842
2843   if (!hinput)
2844   {
2845     AliError(Form("Got a simulation file where I did not find histogram /%s/ALL/EVERYTHING/ANY/INYRANGE/%s",AliAnalysisMuMuBase::MCInputPrefix(),hname.Data()));
2846
2847   }
2848   else
2849   {
2850     r.SetNofInputParticles(*hinput);
2851   }
2852 }
2853
2854 //_____________________________________________________________________________
2855 AliAnalysisMuMuSpectra* AliAnalysisMuMu::SPECTRA(const char* fullpath) const
2856 {
2857   /// Shortcut method to get to a spectra
2858   if (!OC()) return 0x0;
2859
2860   return static_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(fullpath));
2861 }
2862
2863 //_____________________________________________________________________________
2864 void AliAnalysisMuMu::SelectRunByTrigger(const char* triggerList)
2865 {
2866   if (!fMergeableCollection || !fCounterCollection) return;
2867
2868   TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
2869   TIter nextRun(runs);
2870
2871   TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
2872   TIter nextTrigger(triggers);
2873
2874   TObjString* srun;
2875   TObjString* strigger;
2876
2877   TString striggerList(triggerList);
2878
2879   TList* runList = new TList();
2880
2881   while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
2882   {
2883
2884     nextTrigger.Reset();
2885
2886     while ( ( strigger = static_cast<TObjString*>(nextTrigger()) ) )
2887     {
2888       if ( !striggerList.Contains(strigger->String().Data()) )
2889       {
2890         continue;
2891       }
2892
2893       ULong64_t n = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
2894                                                                 strigger->String().Data(),"PSALLHASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00",srun->String().Atoi())));
2895       if ( n > 0 ) runList->Add(srun);
2896
2897     }
2898   }
2899     runList->Sort();
2900     TIter nextRunOK(runList);
2901     while ( ( srun = static_cast<TObjString*>(nextRunOK()) ) )
2902     {
2903       std::cout << srun->String().Atoi() << std::endl;
2904     }
2905     delete runList;
2906
2907     delete triggers;
2908     delete runs;
2909 }
2910 //_____________________________________________________________________________
2911 void AliAnalysisMuMu::TriggerCountCoverage(const char* triggerList,
2912                                            Bool_t compact,
2913                                            Bool_t orderByTriggerCount)
2914 {
2915   // Give the fraction of triggers (in triggerList) relative
2916   // to what is expected in the scalers
2917
2918   TGrid::Connect("alien://"); // to insure the "Trying to connect to server... message does not pollute our output later on...
2919
2920   AliLog::EType_t oldLevel = static_cast<AliLog::EType_t>(AliLog::GetGlobalLogLevel());
2921
2922   AliLog::SetGlobalLogLevel(AliLog::kFatal);
2923
2924   if (!fMergeableCollection || !fCounterCollection) return;
2925
2926   TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
2927   TIter nextRun(runs);
2928
2929   TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
2930   TIter nextTrigger(triggers);
2931
2932   TObjString* srun;
2933   TObjString* strigger;
2934
2935   TString striggerList(triggerList);
2936
2937   ULong64_t total(0);
2938   ULong64_t totalExpected(0);
2939   TString msg;
2940   std::multimap<ULong64_t,std::string> messages;
2941
2942   while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
2943   {
2944     msg.Form("RUN %09d ",srun->String().Atoi());
2945
2946     if (!compact)
2947     {
2948         msg += "\n";
2949     }
2950
2951     ULong64_t nmax(0);
2952
2953     nextTrigger.Reset();
2954
2955     while ( ( strigger = static_cast<TObjString*>(nextTrigger()) ) )
2956     {
2957       if ( !striggerList.Contains(strigger->String().Data()) )
2958       {
2959         continue;
2960       }
2961
2962       ULong64_t n = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
2963                                                             strigger->String().Data(),"ALL",srun->String().Atoi())));
2964
2965       ULong64_t expected = GetTriggerScalerCount(strigger->String().Data(),srun->String().Atoi());
2966
2967
2968       nmax = TMath::Max(n,nmax);
2969
2970       total += n;
2971       totalExpected += expected;
2972
2973       msg += TString::Format("%30s %9lld expected %9lld [%s] ",strigger->String().Data(),n,expected,
2974                              (n>expected ? "!" : " "));
2975
2976       if ( expected > 0 ) {
2977         msg += TString::Format("fraction %5.1f %%",n*100.0/expected);
2978       }
2979
2980       if (!compact)
2981       {
2982         msg += "\n";
2983       }
2984     }
2985     if (nmax>0)
2986     {
2987       if (!orderByTriggerCount)
2988       {
2989         std::cout << msg.Data() << std::endl;
2990       }
2991       else
2992       {
2993         messages.insert(std::make_pair(nmax,static_cast<std::string>(msg.Data())));
2994       }
2995     }
2996   }
2997
2998   std::multimap<ULong64_t,std::string>::const_reverse_iterator it;
2999
3000   ULong64_t current(0);
3001   Int_t n(0);
3002
3003   for ( it = messages.rbegin(); it != messages.rend(); ++it )
3004   {
3005     ++n;
3006     current += it->first;
3007     Double_t percent = ( total > 0.0 ? current*100.0/total : 0.0);
3008     std::cout << Form("%10lld",it->first) << " " << it->second << " percentage of total = " << Form("%7.2f %% %3d",percent,n ) << std::endl;
3009   }
3010
3011   std::cout << Form("--- TOTAL %lld expected %lld fraction %5.1f %%",
3012                     total,totalExpected,totalExpected ? total*100.0/totalExpected : 0.0) << std::endl;
3013
3014
3015
3016   AliLog::SetGlobalLogLevel(oldLevel);
3017   delete triggers;
3018   delete runs;
3019 }
3020
3021 //_____________________________________________________________________________
3022 void AliAnalysisMuMu::UnsetCorrectionPerRun()
3023 {
3024     // drop the correction factors
3025   delete fCorrectionPerRun;
3026   fCorrectionPerRun=0x0;
3027 }
3028
3029 //_____________________________________________________________________________
3030 void AliAnalysisMuMu::Update()
3031 {
3032   /// update the current file with memory
3033
3034   if (!CC() || !OC()) return;
3035
3036   ReOpen(fFilename,"UPDATE");
3037
3038   if (OC())
3039   {
3040     OC()->Write("OC",TObject::kSingleKey|TObject::kOverwrite);
3041   }
3042
3043   ReOpen(fFilename,"READ");
3044
3045   GetCollections(fFilename,fMergeableCollection,fCounterCollection,fBinning,fRunNumbers);
3046 }
3047
3048 //_____________________________________________________________________________
3049 Bool_t AliAnalysisMuMu::Upgrade(const char* filename)
3050 {
3051   /// Upgrade a file
3052   AliAnalysisMuMu m(filename);
3053
3054   return m.Upgrade();
3055 }
3056
3057 //_____________________________________________________________________________
3058 Bool_t AliAnalysisMuMu::Upgrade()
3059 {
3060   /// Upgrade the current file
3061   /// - from single list to one key per object, if needed
3062   /// - from histogramCollection to mergeableCollection, if needed
3063
3064
3065   AliWarning("Out of date method");
3066
3067   TFile* f = ReOpen(fFilename,"UPDATE");
3068
3069   TList* list = static_cast<TList*>(f->Get("chist"));
3070
3071   if (list)
3072   {
3073     // really old file where everything was in a single list
3074
3075     AliHistogramCollection* hc = static_cast<AliHistogramCollection*>(list->At(0));
3076     AliCounterCollection* cc = static_cast<AliCounterCollection*>(list->At(1));
3077
3078     AliMergeableCollection* mc = hc->Convert();
3079
3080     f->cd();
3081
3082     mc->Write("MC",TObject::kSingleKey);
3083     cc->Write("CC",TObject::kSingleKey);
3084
3085     f->Delete("chist;*");
3086
3087     f->Write();
3088
3089   }
3090   else
3091   {
3092     AliHistogramCollection* hc = static_cast<AliHistogramCollection*>(f->Get("HC"));
3093
3094     if ( hc )
3095     {
3096       // old file with histogram collection instead of mergeable collection
3097
3098       AliMergeableCollection* mc = hc->Convert();
3099
3100       f->cd();
3101
3102       mc->Write("MC",TObject::kSingleKey);
3103
3104       f->Delete("HC;*");
3105
3106       f->Write();
3107     }
3108   }
3109
3110   delete f;
3111
3112   return kTRUE;
3113 }
3114
3115 //_____________________________________________________________________________
3116 TH2* AliAnalysisMuMu::ComputeSPDCorrection(const char* type, const char* eventSel, const char* triggerSel, Bool_t bkgReject)
3117 {
3118   TString stype(type);
3119   TString evtype(eventSel);
3120   TString trigtype(triggerSel);
3121
3122   TH2* hNch = static_cast<TH2*>(OC()->Histo(Form("/MCINPUT/%s/%s/V0A/NchVsZVertexVsEta",evtype.Data(),
3123                                                  trigtype.Data()))); // Input Nch // //"/INPUT/QASPDZPSALL/NchVSEtaVSZVertMC"
3124   if ( !hNch )
3125   {
3126     AliError("No Nch histo found");
3127     return 0x0;
3128   }
3129   TH2* hNtr = static_cast<TH2*>(OC()->Histo(Form("/%s/%s/V0A/TrackletsVsZVertexVsEta",evtype.Data(),
3130                                                 trigtype.Data()))); // Reco tracklets //  //"/RECO/QASPDZPSALL/MB1/NtrVSEtaVSZVertMC"
3131   if ( !hNtr )
3132   {
3133     AliError("No tracklets histo found");
3134     return 0x0;
3135   }
3136
3137   TH2* hNtrBkg = static_cast<TH2*>(OC()->Histo(Form("/MCINPUT/%s/%s/V0A/NBkgTrackletsVsZVertexVsEta",evtype.Data(),
3138                                                  trigtype.Data()))); // Reco tracklets //  //"/RECO/QASPDZPSALL/MB1/NtrVSEtaVSZVertMC"
3139   if ( !hNtrBkg )
3140   {
3141     AliError("No background tracklets histo found");
3142     return 0x0;
3143   }
3144
3145
3146   TH2D* hSPDCorr = static_cast<TH2D*>(hNtr->Clone("SPDCorr"));
3147   TString title("");\
3148   if ( stype.Contains("oneOverAccEff")) hSPDCorr->SetTitle("SPD 1/AccxEff correction");
3149   else if ( stype.Contains("AccEffOnly")) hSPDCorr->SetTitle("SPD AccxEff correction");
3150   else if ( stype.Contains("statOneOverAccEff")) hSPDCorr->SetTitle("SPD 1/AccxEff correction stat. unc.");
3151
3152   for (Int_t i = 1 ; i < hNch->GetNbinsX() ; i++)
3153   {
3154     for (Int_t j = 1 ; j < hNch->GetNbinsY() ; j++)
3155     {
3156       Int_t n = hNch->GetBin(i,j);
3157       Double_t nch = hNch->GetBinContent(n);
3158       Double_t ntr = hNtr->GetBinContent(n);
3159       Double_t nBkgtr(0.);
3160       if ( bkgReject ) nBkgtr = hNtrBkg->GetBinContent(n);
3161
3162       Double_t corr(0.),corrErr(0.);
3163       if ( nch != 0. )
3164       {
3165         corr = (ntr - nBkgtr)/nch;
3166         corrErr = TMath::Max( 1./nch,TMath::Sqrt( corr*(1.-corr)/nch ) );
3167       }
3168
3169       if ( stype.Contains("oneOverAccEff"))
3170       {
3171         if ( corr > 0. )
3172         {
3173           hSPDCorr->SetBinContent(n,1./corr);
3174           hSPDCorr->SetBinError(n,corrErr/TMath::Power(corr,2.));
3175         }
3176         else
3177         {
3178           hSPDCorr->SetBinContent(n,0.);
3179           hSPDCorr->SetBinError(n,1.);
3180         }
3181
3182       }
3183       else if ( stype.Contains("AccEffOnly"))
3184       {
3185         hSPDCorr->SetBinContent(n,corr);
3186         hSPDCorr->SetBinError(n,corrErr);
3187       }
3188       else if ( stype.Contains("statOneOverAccEff"))
3189       {
3190         if ( corr != 0. )
3191         {
3192           hSPDCorr->SetBinContent(n,(corrErr/TMath::Power(corr,2.)*100)/(1./corr));
3193         }
3194         else
3195         {
3196           hSPDCorr->SetBinContent(n,-1);
3197         }
3198
3199       }
3200     }
3201   }
3202
3203   return hSPDCorr;
3204 }
3205
3206 //_____________________________________________________________________________
3207 void AliAnalysisMuMu::ComputeFnorm()
3208 {
3209   /// Compute the CMUL to CINT ratio(s)
3210
3211   if (!CC()) return;
3212
3213   OC()->Prune("/FNORM");
3214
3215   AliAnalysisMuMuFnorm computer(*(CC()),AliAnalysisMuMuFnorm::kMUL,Config()->OCDBPath(),Config()->CompactGraphs());
3216
3217   computer.ComputeFnorm();
3218
3219   AliMergeableCollection* fnorm = computer.DetachMC();
3220
3221   OC()->Attach(fnorm,"/FNORM/");
3222
3223   Update();
3224 }
3225
3226 //_____________________________________________________________________________
3227 TH1* AliAnalysisMuMu::ComputeDiffFnormFromHistos(const char* what,const char* quantity,const char* flavour,Bool_t printout)
3228 {
3229   /// Compute the CMUL to CINT ratio(s) from the histos stored in the OC()
3230   //FIXME: This is just a patch...
3231
3232   AliAnalysisMuMuBinning* binning = BIN()->Project(what,quantity,flavour);
3233   if ( !binning )
3234   {
3235     AliError(Form("%s-%s-%s binning does not exist",what,quantity,flavour));
3236     return 0x0;
3237   }
3238   TObjArray* dNchdEtas = binning->CreateBinObjArray();
3239
3240   Double_t* binArray = binning->CreateBinArray();
3241
3242   TIter next(dNchdEtas);
3243   AliAnalysisMuMuBinning::Range* r;
3244
3245   Double_t FNorm(0.);
3246   Double_t FNormError(0.);
3247
3248   TH1* hFNorm = new TH1F("hFNorm","'Global' normalization factor vs dN_{ch}/d#eta;dN_{ch}/d#eta;FNorm",dNchdEtas->GetEntries(),binArray);
3249
3250   Int_t bin(0);
3251   while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(next()) ) )
3252   {
3253
3254     TH1* hCMSL = OC()->Histo(Form("/PSALLHASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00/CMSL7-B-NOPF-MUON/V0A/%s",
3255                              Form("EventsIn%s",r->AsString().Data())));
3256     if ( !hCMSL )
3257     {
3258       AliError(Form("No event histo in bin %s found for CMSL7-B-NOPF-MUON",r->AsString().Data()));
3259       delete binning;
3260       delete dNchdEtas;
3261       delete binArray;
3262       delete hFNorm;
3263       return 0x0;
3264     }
3265
3266     TH1* hCMSLandOMUL = OC()->Histo(Form("/PSALLHASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00/CMSL7-B-NOPF-MUON&0MUL/V0A/%s",
3267                                   Form("EventsIn%s",r->AsString().Data())));
3268     if ( !hCMSLandOMUL )
3269     {
3270       AliError(Form("No event histo in bin %s found for CMSL7-B-NOPF-MUON & 0MUL",r->AsString().Data()));
3271       delete binning;
3272       delete dNchdEtas;
3273       delete binArray;
3274       delete hFNorm;
3275       return 0x0;
3276     }
3277
3278     TH1* hCINT = OC()->Histo(Form("/PSALLHASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00/CINT7-B-NOPF-ALLNOTRD/V0A/%s",
3279                                   Form("EventsIn%s",r->AsString().Data())));
3280     if ( !hCINT )
3281     {
3282       AliError(Form("No event histo in bin %s found for CINT7-B-NOPF-ALLNOTRD",r->AsString().Data()));
3283       delete binning;
3284       delete dNchdEtas;
3285       delete binArray;
3286       delete hFNorm;
3287       return 0x0;
3288     }
3289
3290     TH1* hCINTandOMSL = OC()->Histo(Form("/PSALLHASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00/CINT7-B-NOPF-ALLNOTRD&0MSL/V0A/%s",
3291                                   Form("EventsIn%s",r->AsString().Data())));
3292     if ( !hCINTandOMSL )
3293     {
3294       AliError(Form("No event histo in bin %s found for CINT7-B-NOPF-ALLNOTRD & 0MSL",r->AsString().Data()));
3295       delete binning;
3296       delete dNchdEtas;
3297       delete binArray;
3298       delete hFNorm;
3299       return 0x0;
3300     }
3301
3302     FNorm = (hCMSL->GetBinContent(1)/hCMSLandOMUL->GetBinContent(1))*(hCINT->GetBinContent(1)/hCINTandOMSL->GetBinContent(1));
3303     FNormError = ErrorPropagationAxBoverCxD(hCMSL->GetBinContent(1),hCINT->GetBinContent(1),hCMSLandOMUL->GetBinContent(1),hCINTandOMSL->GetBinContent(1));
3304
3305     if ( printout ) std::cout << r->AsString().Data() << " : " << FNorm << " +- " << FNormError << std::endl;
3306
3307     hFNorm->SetBinContent(++bin,FNorm);
3308     hFNorm->SetBinError(bin,FNormError);
3309   }
3310
3311   delete binning;
3312   delete dNchdEtas;
3313   delete[] binArray;
3314
3315   return hFNorm;
3316 }
3317
3318 //_____________________________________________________________________________
3319 void AliAnalysisMuMu::ComputeDiffFnormFromInt(const char* triggerCluster, const char* eventSelection, AliMergeableCollection* mc, const char* what,const char* quantity,const char* flavour,Bool_t printout)
3320 {
3321   TString striggerCluster(triggerCluster);
3322   if ( striggerCluster.Contains("MUON") && !striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON";
3323   else if ( striggerCluster.Contains("ALLNOTRD") && !striggerCluster.Contains("MUON") ) striggerCluster = "ALLNOTRD";
3324   else if ( striggerCluster.Contains("MUON") && striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON-ALLNOTRD";
3325   else
3326   {
3327     AliError("Unknown trigger cluster");
3328     return;
3329   }
3330
3331   TString seventSelection(eventSelection);
3332
3333   TString id(Form("/FNORM-%s/%s/V0A",striggerCluster.Data(),seventSelection.Data()));
3334
3335   AliAnalysisMuMuBinning* binning = BIN()->Project(what,quantity,flavour);
3336   if ( !binning )
3337   {
3338     AliError(Form("%s-%s-%s binning does not exist",what,quantity,flavour));
3339     return;
3340   }
3341
3342   TString path(Form("%s/%s/%s",
3343                     First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kFALSE)).Data(),
3344                     First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data(),
3345                     First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,kFALSE)).Data()));
3346   if ( !mc )
3347   {
3348     AliError("Error: No mergeable collection to get Nch histo");
3349     delete binning;
3350     return;
3351   }
3352   TH1* hNtr = static_cast<TH1*>(mc->Histo(Form("/%s/Nch",path.Data())));
3353   if ( !hNtr )
3354   {
3355     AliError(Form("Error: No /%s/Nch histo in mergeable collection",path.Data()));
3356     delete binning;
3357     return;
3358   }
3359   Int_t nTrackletsCorrTot = hNtr->Integral();
3360
3361   TObjArray* bin = binning->CreateBinObjArray(what,quantity,flavour);
3362   Int_t nEntries = bin->GetEntries();
3363   Double_t* binArray = binning->CreateBinArray();
3364   Double_t FNormTot(0.);
3365   Double_t FNormTotError(0.);
3366
3367   TH1* hNorm = OC()->Histo(Form("%s/hFNormInt",id.Data()));
3368
3369   TH1* hFNormTot = new TH1F("hFNormVSdNchdEtaFromInt","Normalization factor vs dN_{ch}/d#eta;dN_{ch}/d#eta;FNorm",nEntries,binArray);
3370
3371   Double_t FNormGlobal = hNorm->GetBinContent(1);
3372   Double_t FNormGlobalError = hNorm->GetBinError(1);
3373
3374   if ( printout ) std::cout << "Global FNorm = " << FNormGlobal << " + - " << FNormGlobalError << std::endl;
3375
3376   TIter nextBin(bin);
3377   AliAnalysisMuMuBinning::Range* r;
3378   Int_t i(1);
3379   while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) ) //Bin loop
3380   {
3381     Double_t nTrklsCorrBin = hNtr->Integral(r->Xmin(),r->Xmax());
3382     Double_t nTrklsCorrBinFrac = nTrklsCorrBin / nTrackletsCorrTot;
3383     Double_t nTrklsCorrBinFracError = TMath::Sqrt( TMath::Power(TMath::Sqrt(nTrklsCorrBin) / nTrackletsCorrTot,2.) +
3384                                                   TMath::Power(TMath::Sqrt(nTrackletsCorrTot)*nTrklsCorrBin / TMath::Power(nTrackletsCorrTot,2.) ,2.) );
3385
3386     FNormTot = FNormGlobal*nTrklsCorrBinFrac;
3387     FNormTotError = TMath::Sqrt( TMath::Power(FNormGlobalError*nTrklsCorrBinFrac,2.) + TMath::Power(FNormGlobal*nTrklsCorrBinFracError,2.) );
3388
3389     hFNormTot->SetBinContent(i,FNormTot);
3390     hFNormTot->SetBinError(i,FNormTotError);
3391     i++;
3392
3393     if ( printout ) std::cout << "Bin: " << r->AsString().Data() << " ; " << " FNorm = " << FNormTot << " +- " << FNormTotError << std::endl;
3394   }
3395
3396   TH1* o = fMergeableCollection->Histo(id.Data(),hFNormTot->GetName());
3397
3398   if (o)
3399   {
3400     AliWarning(Form("Replacing %s/%s",id.Data(),hFNormTot->GetName()));
3401     fMergeableCollection->Remove(Form("%s/%s",id.Data(),hFNormTot->GetName()));
3402   }
3403
3404   Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),hFNormTot);
3405
3406   if ( adoptOK ) std::cout << "+++FNorm histo " << hFNormTot->GetName() << " adopted" << std::endl;
3407   else AliError(Form("Could not adopt FNorm histo %s",hFNormTot->GetName()));
3408
3409   delete binning;
3410   delete bin;
3411   delete[] binArray;
3412 }
3413
3414 //_____________________________________________________________________________
3415 void AliAnalysisMuMu::ComputeDiffFnormFromCounters(const char* triggerCluster, const char* eventSelection, const char* filePileUpCorr, const char* what,const char* quantity,const char* flavour,Bool_t printout)
3416 {
3417   /// Compute the CMUL to CINT ratio(s) in 2 steps from the CC(), in bins
3418
3419   TString colType(First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data());
3420   if ( colType.Contains("-B-") ) colType = "B";
3421   else if ( colType.Contains("-S-") ) colType = "S";
3422   else
3423   {
3424     AliError("Unknown collision type");
3425     return;
3426   }
3427
3428   TString triggerType(First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data());
3429   if ( triggerType.Contains("7-") ) triggerType = "7";
3430   else if ( triggerType.Contains("8-") ) triggerType = "8";
3431   else
3432   {
3433     AliError("Unknown trigger type");
3434     return;
3435   }
3436
3437   TString striggerCluster(triggerCluster);
3438   if ( striggerCluster.Contains("MUON") && !striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON";
3439   else if ( striggerCluster.Contains("ALLNOTRD") && !striggerCluster.Contains("MUON") ) striggerCluster = "ALLNOTRD";
3440   else if ( striggerCluster.Contains("MUON") && striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON-ALLNOTRD";
3441   else
3442   {
3443     AliError("Unknown trigger cluster");
3444     return;
3445   }
3446
3447   std::cout << striggerCluster.Data() << std::endl;
3448
3449   Bool_t corrPU(kFALSE);
3450   TObjArray* pUCorr = new TObjArray();
3451   if ( strlen(filePileUpCorr) > 0 )
3452   {
3453     //    std::cout << "Extracting Pile-Up correction factors from " << filePileUpCorr << std::endl;
3454     char line[1024];
3455     ifstream in(filePileUpCorr);
3456
3457     while ( in.getline(line,1024,'\n'))
3458     {
3459       TString lrun(line);
3460       TString lvalue(line);
3461
3462       lrun.Remove(0,4);
3463       lrun.Remove(6,67);
3464
3465       lvalue.Remove(0,57);//71
3466
3467       //      std::cout << "RUN: " << lrun.Data() << " PUFactor = " << lvalue.Data() << std::endl;
3468
3469       pUCorr->Add(new TParameter<Double_t>(lrun.Data(),lvalue.Atof()));
3470     }
3471     corrPU = kTRUE;
3472   }
3473
3474   TString seventSelection(eventSelection);
3475   TString sruns = CC()->GetKeyWords("run");
3476   TObjArray* runs = sruns.Tokenize(",");
3477   Double_t NofRuns = runs->GetEntries();
3478
3479   TIter nextRun(runs);
3480   TObjString* s;
3481
3482   AliAnalysisMuMuBinning* binning = BIN()->Project(what,quantity,flavour);
3483   if ( !binning )
3484   {
3485     AliError(Form("%s-%s-%s binning does not exist",what,quantity,flavour));
3486     return;
3487   }
3488   TObjArray* bin = binning->CreateBinObjArray(what,quantity,flavour);
3489   Double_t* binArray = binning->CreateBinArray();
3490   Int_t nEntries = bin->GetEntries();
3491
3492   TH1* h;
3493   TH1* hNofEqMB = new TH1F("hNofEqMBVSdNchdEta","Equivalent MB events per CMUL for vs dN_{ch}/d#eta",bin->GetEntries(),binArray);
3494   TH1* hFNormTot = new TH1F("hFNormVSdNchdEta","Normalization factor vs dN_{ch}/d#eta;dN_{ch}/d#eta;FNorm",bin->GetEntries(),binArray);
3495
3496   Double_t* FNormTot = new Double_t[nEntries];
3497   Double_t* FNormTotError = new Double_t[nEntries];
3498
3499   TString id(Form("/FNORM-%s/%s/V0A",striggerCluster.Data(),seventSelection.Data()));//HASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00
3500
3501   TList* lRun2Reject = new TList();
3502   lRun2Reject->SetOwner(kTRUE);
3503
3504   Int_t i(0); //dNchdEta bin number
3505   TObjArray* aCluster = striggerCluster.Tokenize("-");
3506   TIter nextCluster(aCluster);
3507   TObjString* striggerClusterS;
3508
3509   TIter nextBin(bin);
3510   AliAnalysisMuMuBinning::Range* r;
3511   while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) ) //Bin loop
3512   {
3513     FNormTot[i] = 0;
3514     FNormTotError[i] = 0;
3515     if ( printout )
3516     {
3517       std::cout << "______________________________" << std::endl;
3518       std::cout << "Bin: " << r->AsString().Data() << std::endl;
3519     }
3520
3521     h = new TH1F(Form("hFNormVSrun_%s",r->AsString().Data()),Form("Normalization factor vs run for %s ;run;FNorm",r->AsString().Data()),NofRuns,1,NofRuns);
3522     //Set the run labels
3523
3524     Double_t nCMULBin = CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A/bin:%s",
3525                                           seventSelection.Data(),triggerType.Data(),colType.Data(),r->AsString().Data())); //Nof CMUL7/8 events in Bin summed over runs
3526    //HASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00
3527
3528     Int_t j(1); //Run label index
3529     nextRun.Reset();
3530     while ( ( s = static_cast<TObjString*>(nextRun())) ) //Run loop
3531     {
3532       Double_t nCMSL(0.),nCMSLandOMUL(0.);
3533       nextCluster.Reset();
3534       while ( (striggerClusterS = static_cast<TObjString*>(nextCluster())) && nCMSL == 0. )
3535       {
3536         nCMSL = CC()->GetSum(Form("/event:%s/trigger:CMSL%s-%s-NOPF-%s/centrality:V0A/run:%s/bin:%s",
3537                                   seventSelection.Data(),triggerType.Data(),colType.Data(),striggerClusterS->GetName(),s->GetName(),r->AsString().Data()));
3538
3539         nCMSLandOMUL = CC()->GetSum(Form("/event:%s/trigger:CMSL%s-%s-NOPF-%s&0MUL/centrality:V0A/run:%s/bin:%s",
3540                                          seventSelection.Data(),triggerType.Data(),colType.Data(),striggerClusterS->GetName(),s->GetName(),r->AsString().Data()));
3541       }
3542       Double_t nCMUL = CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A/run:%s/bin:%s",
3543                                 seventSelection.Data(),triggerType.Data(),colType.Data(),s->GetName(),r->AsString().Data()));
3544
3545       Double_t nCINT = CC()->GetSum(Form("/event:%s/trigger:CINT%s-%s-NOPF-ALLNOTRD/centrality:V0A/run:%s/bin:%s",
3546                                          seventSelection.Data(),triggerType.Data(),colType.Data(),s->GetName(),r->AsString().Data()));
3547
3548       Double_t nCINTandOMSL = CC()->GetSum(Form("/event:%s/trigger:CINT%s-%s-NOPF-ALLNOTRD&0MSL/centrality:V0A/run:%s/bin:%s",
3549                                                 seventSelection.Data(),triggerType.Data(),colType.Data(),s->GetName(),r->AsString().Data()));
3550
3551       Double_t FNorm(0.);
3552       Double_t FNormError(0.);
3553       Double_t FNormError2(0.);
3554       Double_t pUfactor = 1.;
3555       if ( nCMSLandOMUL != 0. && nCINTandOMSL !=0. && nCMSL != 0. && nCINT !=0. )
3556       {
3557         if (corrPU)
3558         {
3559           TParameter<Double_t>* p = static_cast<TParameter<Double_t>*>(pUCorr->FindObject(s->GetName()));
3560           if ( p ) pUfactor = p->GetVal();
3561           else
3562           {
3563             AliError(Form("Run %s not found in pile-up correction list",s->GetName()));
3564           }
3565         }
3566
3567         FNorm = (nCMSL*nCINT)*pUfactor/(nCMSLandOMUL*nCINTandOMSL);
3568         FNormError = ErrorPropagationAxBoverCxD(nCMSL,nCINT,nCMSLandOMUL,nCINTandOMSL)*pUfactor;
3569         FNormError2 = AliAnalysisMuMuResult::ErrorABCD(nCMSL, TMath::Sqrt(nCMSL), nCINT, TMath::Sqrt(nCINT), nCMSLandOMUL, TMath::Sqrt(nCMSLandOMUL), nCINTandOMSL, TMath::Sqrt(nCINTandOMSL));
3570       }
3571       else
3572       {
3573         if ( nCINT == 0 ) std::cout << " Warning: Run " << s->GetName() << " has no MB trigger in this bin" << std::endl;
3574
3575         lRun2Reject->Add(new TObjString(s->GetName()));
3576         if ( printout ) std::cout << "Run " << s->GetName() << " not used for FNorm cause lack of stats" << std::endl;
3577         continue;
3578       }
3579       FNormTot[i] += FNorm*nCMUL; // This is the sum of equivalent Nof MB per CMUL run by run. NOTE: This sum is NOT always the total equivalent Nof MB per CMUL because in pp 2012 if just one cluster is used at a time this sum is not the sum for all runs
3580       FNormTotError[i] += TMath::Power(nCMUL*FNormError,2.) + TMath::Power(FNorm*TMath::Sqrt(nCMUL),2.);
3581
3582       if ( printout ) std::cout << "Run " << s->GetName() << " FNorm = " << FNorm << " +- " << FNormError << " (" << FNormError2 << ")" << " ; PUFactor =" << pUfactor << " ; " << "Nof CMUL = " << nCMUL << std::endl;
3583
3584       h->GetXaxis()->SetBinLabel(j,s->GetName());
3585       h->SetBinContent(j,FNorm);
3586       h->SetBinError(j++,FNormError);
3587
3588     }
3589
3590 //    std::cout << "NofCMUL in " << i << " = " << nCMULBin << std::endl;
3591
3592     TIter nextRejectRun(lRun2Reject);
3593     TObjString* run2Rej;
3594     Double_t nCMULBinRej(0.);
3595     while ( (run2Rej = static_cast<TObjString*>(nextRejectRun())) )
3596     {
3597       nCMULBinRej += CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A/bin:%s/run:%s",
3598                                        seventSelection.Data(),triggerType.Data(),colType.Data(),r->AsString().Data(),
3599                                        run2Rej->GetName())); //Sum of CMUL7 events from rejected runs
3600     }
3601
3602     nCMULBin = nCMULBin - nCMULBinRej;
3603     lRun2Reject->Clear();
3604
3605     FNormTotError[i] =  TMath::Sqrt(TMath::Power(TMath::Sqrt(FNormTotError[i])/nCMULBin,2.) + TMath::Power(FNormTot[i]*TMath::Sqrt(nCMULBin)/TMath::Power(nCMULBin,2.),2.));
3606     FNormTot[i] = FNormTot[i]/nCMULBin;
3607
3608     std::cout << "Mean FNorm in Bin = " << FNormTot[i]  << " +- " << FNormTotError[i] <<std::endl;
3609
3610     hFNormTot->SetBinContent(i+1,FNormTot[i]);
3611     hFNormTot->SetBinError(i+1,FNormTotError[i]);
3612
3613     //____
3614     nCMULBin = CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A/bin:%s",
3615                                  seventSelection.Data(),triggerType.Data(),colType.Data(),r->AsString().Data())); //Nof CMUL7/8 events in Bin summed over runs
3616
3617     Double_t nofEqMB = FNormTot[i]*nCMULBin;
3618     Double_t nofEqMBError = TMath::Sqrt( TMath::Power(FNormTotError[i]*nCMULBin,2.) + TMath::Power(FNormTot[i]*TMath::Sqrt(nCMULBin),2.) );
3619
3620     std::cout << "EqMB in Bin  = " << nofEqMB << " +- " << nofEqMBError << " ; nCMUL = " << nCMULBin << std::endl;
3621
3622     hNofEqMB->SetBinContent(i+1,nofEqMB);
3623     hNofEqMB->SetBinError(i+1,nofEqMBError);
3624     //____
3625
3626     TH1* o = fMergeableCollection->Histo(id.Data(),h->GetName());
3627
3628     if (o)
3629     {
3630       AliWarning(Form("Replacing %s/%s",id.Data(),h->GetName()));
3631       fMergeableCollection->Remove(Form("%s/%s",id.Data(),h->GetName()));
3632     }
3633
3634     Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),h);
3635
3636     if ( adoptOK ) std::cout << "+++FNorm histo " << h->GetName() << " adopted" << std::endl;
3637     else AliError(Form("Could not adopt FNorm histo %s",h->GetName()));
3638
3639     i++;
3640     lRun2Reject->Clear();
3641   }
3642
3643   TH1* o = fMergeableCollection->Histo(id.Data(),hFNormTot->GetName());
3644
3645   if (o)
3646   {
3647     AliWarning(Form("Replacing %s/%s",id.Data(),hFNormTot->GetName()));
3648     fMergeableCollection->Remove(Form("%s/%s",id.Data(),hFNormTot->GetName()));
3649   }
3650
3651   Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),hFNormTot);
3652
3653   if ( adoptOK ) std::cout << "+++FNorm histo " << hFNormTot->GetName() << " adopted" << std::endl;
3654   else AliError(Form("Could not adopt FNorm histo %s",hFNormTot->GetName()));
3655
3656
3657   o = fMergeableCollection->Histo(id.Data(),hNofEqMB->GetName());
3658
3659   if (o)
3660   {
3661     AliWarning(Form("Replacing %s/%s",id.Data(),hNofEqMB->GetName()));
3662     fMergeableCollection->Remove(Form("%s/%s",id.Data(),hNofEqMB->GetName()));
3663   }
3664
3665   adoptOK = fMergeableCollection->Adopt(id.Data(),hNofEqMB);
3666
3667   if ( adoptOK ) std::cout << "+++FNorm histo " << hNofEqMB->GetName() << " adopted" << std::endl;
3668   else AliError(Form("Could not adopt FNorm histo %s",hNofEqMB->GetName()));
3669
3670   delete binning;
3671   delete runs;
3672   delete aCluster;
3673   delete bin;
3674   delete[] binArray;
3675   delete[] FNormTot;
3676   delete[] FNormTotError;
3677   delete lRun2Reject;
3678
3679   return;
3680
3681 }
3682
3683 //_____________________________________________________________________________
3684 void AliAnalysisMuMu::ComputeDiffFnormFromGlobal(const char* triggerCluster, const char* eventSelection, const char* what,const char* quantity,
3685                                                  const char* flavour, Bool_t printout)
3686 {
3687   TString colType(First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data());
3688   if ( colType.Contains("-B-") ) colType = "B";
3689   else if ( colType.Contains("-S-") ) colType = "S";
3690   else
3691   {
3692     AliError("Unknown collision type");
3693     return;
3694   }
3695
3696   TString triggerType(First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data());
3697   if ( triggerType.Contains("7-") ) triggerType = "7";
3698   else if ( triggerType.Contains("8-") ) triggerType = "8";
3699   else
3700   {
3701     AliError("Unknown trigger type");
3702     return;
3703   }
3704
3705   TString striggerCluster(triggerCluster);
3706   if ( striggerCluster.Contains("MUON") && !striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON";
3707   else if ( striggerCluster.Contains("ALLNOTRD") && !striggerCluster.Contains("MUON") ) striggerCluster = "ALLNOTRD";
3708   else if ( striggerCluster.Contains("MUON") && striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON-ALLNOTRD";
3709   else
3710   {
3711     AliError("Unknown trigger cluster");
3712     return;
3713   }
3714
3715   TString seventSelection(eventSelection);
3716   TString id(Form("/FNORM-%s/%s/V0A",striggerCluster.Data(),seventSelection.Data()));
3717
3718   TH1* hFnormGlobal = OC()->Histo(id.Data(),"hFNormVSdNchdEta");
3719   if( !hFnormGlobal)
3720   {
3721     AliError("hFNormVSdNchdEta not found");
3722     return;
3723   }
3724
3725   TH1* hFnormGlobalInt = OC()->Histo(id.Data(),"hFNormInt");
3726   if( !hFnormGlobalInt)
3727   {
3728     AliError("hFNormInt not found");
3729     return;
3730   }
3731   Double_t FNormGlobal = hFnormGlobalInt->GetBinContent(1);
3732   Double_t FNormGlobalError = hFnormGlobalInt->GetBinError(1);
3733
3734   AliAnalysisMuMuBinning* binning = BIN()->Project(what,quantity,flavour);
3735   if ( !binning )
3736   {
3737     AliError(Form("%s-%s-%s binning does not exist",what,quantity,flavour));
3738     return;
3739   }
3740   TObjArray* bin = binning->CreateBinObjArray(what,quantity,flavour);
3741
3742   TH1* hFNormVSNtr = static_cast<TH1*>(hFnormGlobal->Clone());
3743   hFNormVSNtr->SetTitle("Normalization factor vs dN_{ch}/d#eta;dN_{ch}/d#eta;FNorm");
3744   hFNormVSNtr->SetName("hFNormVSdNchdEtaFromGlobal");
3745
3746   TH1* hNMBVSNtr = static_cast<TH1*>(hFnormGlobal->Clone());
3747   hNMBVSNtr->SetTitle("Equivalent MB events per CMUL for vs dN_{ch}/d#eta;NMB");
3748   hNMBVSNtr->SetName("hNofEqMBVSdNchdEtaFromGlobal");
3749
3750
3751   Double_t nCMULTot = CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A",
3752                                         seventSelection.Data(),triggerType.Data(),colType.Data()));
3753
3754   Double_t nCINTTot = CC()->GetSum(Form("/event:%s/trigger:CINT%s-%s-NOPF-ALLNOTRD/centrality:V0A",
3755                                      seventSelection.Data(),triggerType.Data(),colType.Data()));
3756   TIter nextBin(bin);
3757   AliAnalysisMuMuBinning::Range* r;
3758   Int_t i(1);
3759   while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) ) //Bin loop
3760   {
3761     Double_t nCMUL = CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A/bin:%s",
3762                                        seventSelection.Data(),triggerType.Data(),colType.Data(),r->AsString().Data()));
3763
3764     Double_t nCINT = CC()->GetSum(Form("/event:%s/trigger:CINT%s-%s-NOPF-ALLNOTRD/centrality:V0A/bin:%s",
3765                                        seventSelection.Data(),triggerType.Data(),colType.Data(),r->AsString().Data()));
3766
3767     Double_t f = nCMUL/nCMULTot;
3768     Double_t fError = TMath::Sqrt( TMath::Power(TMath::Sqrt(nCMUL)/nCMULTot,2.) + TMath::Power(nCMUL*TMath::Sqrt(nCMULTot)/TMath::Power(nCMULTot,2.),2.) );
3769
3770     Double_t g = nCINT/nCINTTot;
3771     Double_t gError = TMath::Sqrt( TMath::Power(TMath::Sqrt(nCINT)/nCINTTot,2.) + TMath::Power(nCINT*TMath::Sqrt(nCINTTot)/TMath::Power(nCINTTot,2.),2.) );
3772
3773     Double_t value = FNormGlobal*(g/f);
3774     Double_t error = TMath::Sqrt( TMath::Power(FNormGlobalError*(g/f),2.) + TMath::Power(FNormGlobal*(gError/f),2.) + TMath::Power(FNormGlobal*g*fError/TMath::Power(f,2.),2.) );
3775
3776     hFNormVSNtr->SetBinContent(i,value);
3777     hFNormVSNtr->SetBinError(i,error);
3778
3779     if (printout) std::cout << value << " +- " << error << " ; nCMUL = " << nCMUL << std::endl;
3780
3781     hNMBVSNtr->SetBinContent(i,value*nCMUL);
3782     hNMBVSNtr->SetBinError(i,TMath::Sqrt( TMath::Power(error*nCMUL,2.) + TMath::Power(value*TMath::Sqrt(nCMUL),2.) ));
3783
3784     i++;
3785   }
3786
3787   TH1* o = fMergeableCollection->Histo(id.Data(),hFNormVSNtr->GetName());
3788
3789   if (o)
3790   {
3791     AliWarning(Form("Replacing %s/%s",id.Data(),hFNormVSNtr->GetName()));
3792     fMergeableCollection->Remove(Form("%s/%s",id.Data(),hFNormVSNtr->GetName()));
3793   }
3794
3795   Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),hFNormVSNtr);
3796
3797   if ( adoptOK ) std::cout << "+++FNorm histo " << hFNormVSNtr->GetName() << " adopted" << std::endl;
3798   else AliError(Form("Could not adopt FNorm histo %s",hFNormVSNtr->GetName()));
3799
3800   o = fMergeableCollection->Histo(id.Data(),hNMBVSNtr->GetName());
3801
3802   if (o)
3803   {
3804     AliWarning(Form("Replacing %s/%s",id.Data(),hNMBVSNtr->GetName()));
3805     fMergeableCollection->Remove(Form("%s/%s",id.Data(),hNMBVSNtr->GetName()));
3806   }
3807
3808   adoptOK = fMergeableCollection->Adopt(id.Data(),hNMBVSNtr);
3809
3810   if ( adoptOK ) std::cout << "+++FNorm histo " << hNMBVSNtr->GetName() << " adopted" << std::endl;
3811   else AliError(Form("Could not adopt FNorm histo %s",hNMBVSNtr->GetName()));
3812
3813
3814 }
3815
3816 //_____________________________________________________________________________
3817 void AliAnalysisMuMu::ComputeMeanFnorm(const char* triggerCluster, const char* eventSelection, const char* what,const char* quantity,const char* flavour, Bool_t printout)
3818 {
3819   /// Compute the mean Fnorm and mean NMB from the offline and "rescaled global" methods
3820
3821   TString seventSelection(eventSelection);
3822   TString striggerCluster(triggerCluster);
3823   if ( striggerCluster.Contains("MUON") && !striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON";
3824   else if ( striggerCluster.Contains("ALLNOTRD") && !striggerCluster.Contains("MUON") ) striggerCluster = "ALLNOTRD";
3825   else if ( striggerCluster.Contains("MUON") && striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON-ALLNOTRD";
3826   else
3827   {
3828     AliError("Unknown trigger cluster");
3829     return;
3830   }
3831
3832   TString triggerType(First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data());
3833   if ( triggerType.Contains("7-") ) triggerType = "7";
3834   else if ( triggerType.Contains("8-") ) triggerType = "8";
3835   else
3836   {
3837     AliError("Unknown trigger type");
3838     return;
3839   }
3840
3841   TString colType(First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data());
3842   if ( colType.Contains("-B-") ) colType = "B";
3843   else if ( colType.Contains("-S-") ) colType = "S";
3844   else
3845   {
3846     AliError("Unknown collision type");
3847     return;
3848   }
3849
3850   TH1* hMB = OC()->Histo(Form("/FNORM-%s/%s/V0A/hNofEqMBVSdNchdEta",striggerCluster.Data(),seventSelection.Data()));//HASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00
3851   if ( !hMB )
3852   {
3853     AliError("Histo hNofEqMBVSdNchdEta not found");
3854     return;
3855   }
3856
3857   TH1* hMBG = OC()->Histo(Form("/FNORM-%s/%s/V0A/hNofEqMBVSdNchdEtaFromGlobal",striggerCluster.Data(),seventSelection.Data()));//HASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00
3858   if ( !hMBG )
3859   {
3860     AliError("Histo hNofEqMBVSdNchdEtaFromGlobal not found");
3861     return;
3862   }
3863
3864   AliAnalysisMuMuBinning* binning = BIN()->Project(what,quantity,flavour);
3865   if ( !binning )
3866   {
3867     AliError(Form("%s-%s-%s binning does not exist",what,quantity,flavour));
3868     return;
3869   }
3870   TObjArray* bin = binning->CreateBinObjArray(what,quantity,flavour);
3871
3872
3873   TString id(Form("/FNORM-%s/%s/V0A",striggerCluster.Data(),seventSelection.Data()));
3874
3875   TH1* hMBMean = static_cast<TH1*>(hMBG->Clone());
3876   hMBMean->SetName("hNofEqMBVSdNchdEtaFromMean");
3877
3878   TH1* hFnormMean = static_cast<TH1*>(hMBG->Clone());
3879   hFnormMean->SetName("hFNormVSdNchdEtaFromMean");
3880
3881   for ( Int_t i = 1 ; i <= hMB->GetNbinsX() ; i++ )
3882   {
3883     Double_t Fn = hMB->GetBinContent(i);
3884     Double_t Fng = hMBG->GetBinContent(i);
3885
3886     Double_t FnE = hMB->GetBinError(i);
3887     Double_t FngE = hMBG->GetBinError(i);
3888
3889     Double_t meanBin = (Fn + Fng) / 2.;
3890     Double_t meanBinError = TMath::Sqrt( TMath::Power(FnE/2.,2.) + TMath::Power(FngE/2.,2.) );
3891     Double_t meanBinSys = TMath::Abs( meanBin - Fn );
3892
3893 //    Double_t meanBin = (Fn/TMath::Power(FnE,2.) + Fng/TMath::Power(FngE,2.)) / ( 1./TMath::Power(FnE,2.) + 1./TMath::Power(FngE,2.) );
3894 //    Double_t meanBinError = 1. / TMath::Sqrt( 1./TMath::Power(FnE,2.) + 1./TMath::Power(FngE,2.) );
3895 //    Double_t meanBinSys = TMath::Sqrt( TMath::Power(Fn - meanBin,2.)/TMath::Power(FnE,2.) + TMath::Power(Fng - meanBin,2.)/TMath::Power(FngE,2.) );
3896
3897
3898     hMBMean->SetBinContent(i,meanBin);
3899     hMBMean->SetBinError(i,meanBinError);
3900
3901     Double_t nCMULBin = CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A/bin:%s",
3902                                           seventSelection.Data(),triggerType.Data(),colType.Data(),
3903                                           static_cast<AliAnalysisMuMuBinning::Range*>(bin->At(i-1))->AsString().Data()));
3904
3905     if (printout) std::cout << meanBinSys/nCMULBin << std::endl;
3906
3907     Double_t meanFnBin = meanBin/nCMULBin;
3908     Double_t meanFnBinError = TMath::Sqrt( TMath::Power(meanBinError/nCMULBin,2) + TMath::Power(meanBin/TMath::Power(nCMULBin,2.),2) );
3909
3910     if (printout) std::cout << meanBinSys/nCMULBin/meanFnBin << std::endl;
3911
3912     if (printout) std::cout << meanFnBin << " +- " << meanFnBinError << std::endl;
3913
3914     hFnormMean->SetBinContent(i,meanFnBin);
3915     hFnormMean->SetBinError(i,meanFnBinError);
3916   }
3917
3918   TH1* o = fMergeableCollection->Histo(id.Data(),hMBMean->GetName());
3919
3920   if (o)
3921   {
3922     AliWarning(Form("Replacing %s/%s",id.Data(),hMBMean->GetName()));
3923     fMergeableCollection->Remove(Form("%s/%s",id.Data(),hMBMean->GetName()));
3924   }
3925
3926   Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),hMBMean);
3927
3928   if ( adoptOK ) std::cout << "+++FNorm histo " << hMBMean->GetName() << " adopted" << std::endl;
3929   else AliError(Form("Could not adopt FNorm histo %s",hMBMean->GetName()));
3930
3931
3932   o = fMergeableCollection->Histo(id.Data(),hFnormMean->GetName());
3933
3934   if (o)
3935   {
3936     AliWarning(Form("Replacing %s/%s",id.Data(),hFnormMean->GetName()));
3937     fMergeableCollection->Remove(Form("%s/%s",id.Data(),hFnormMean->GetName()));
3938   }
3939
3940   adoptOK = fMergeableCollection->Adopt(id.Data(),hFnormMean);
3941
3942   if ( adoptOK ) std::cout << "+++FNorm histo " << hFnormMean->GetName() << " adopted" << std::endl;
3943   else AliError(Form("Could not adopt FNorm histo %s",hFnormMean->GetName()));
3944
3945
3946 }
3947
3948 //_____________________________________________________________________________
3949 void AliAnalysisMuMu::ComputeIntFnormFromCounters(const char* triggerCluster, const char* eventSelection, const char* filePileUpCorr, Bool_t printout)
3950 {
3951   /// Compute the CMUL to CINT ratio(s) in 2 steps from the CC(), integrated
3952
3953   TString colType(First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data());
3954   if ( colType.Contains("-B-") ) colType = "B";
3955   else if ( colType.Contains("-S-") ) colType = "S";
3956   else
3957   {
3958     AliError("Unknown collision type");
3959     return;
3960   }
3961
3962   TString triggerType(First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data());
3963   if ( triggerType.Contains("7-") ) triggerType = "7";
3964   else if ( triggerType.Contains("8-") ) triggerType = "8";
3965   else
3966   {
3967     AliError("Unknown trigger type");
3968     return;
3969   }
3970
3971   TString striggerCluster(triggerCluster);
3972   if ( striggerCluster.Contains("MUON") && !striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON";
3973   else if ( striggerCluster.Contains("ALLNOTRD") && !striggerCluster.Contains("MUON") ) striggerCluster = "ALLNOTRD";
3974   else if ( striggerCluster.Contains("MUON") && striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON-ALLNOTRD";
3975   else
3976   {
3977     AliError("Unknown trigger cluster");
3978     return;
3979   }
3980
3981   TString seventSelection(eventSelection);
3982   TString sruns = CC()->GetKeyWords("run");
3983   TObjArray* runs = sruns.Tokenize(",");
3984   Double_t NofRuns = runs->GetEntries();
3985
3986   Bool_t corrPU(kFALSE);
3987   TObjArray* pUCorr = new TObjArray();
3988   if ( strlen(filePileUpCorr) > 0 )
3989   {
3990 //    std::cout << "Extracting Pile-Up correction factors from " << filePileUpCorr << std::endl;
3991     char line[1024];
3992     ifstream in(filePileUpCorr);
3993
3994     while ( in.getline(line,1024,'\n'))
3995     {
3996       TString lrun(line);
3997       TString lvalue(line);
3998
3999       lrun.Remove(0,4);
4000       lrun.Remove(6,67);
4001
4002       lvalue.Remove(0,57);//71
4003
4004 //      std::cout << "RUN: " << lrun.Data() << " PUFactor = " << lvalue.Data() << std::endl;
4005
4006       pUCorr->Add(new TParameter<Double_t>(lrun.Data(),lvalue.Atof()));
4007     }
4008     corrPU = kTRUE;
4009   }
4010
4011   TIter nextRun(runs);
4012   TObjString* s;
4013
4014   TH1* h;
4015
4016   Double_t FNormTot(0.);
4017   Double_t FNormTotError(0.);
4018
4019   TString id(Form("/FNORM-%s/%s/V0A",striggerCluster.Data(),seventSelection.Data()));//HASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00
4020
4021   h = new TH1F("hFNormIntVSrun","Integrated Normalization factor vs run;run;FNorm",NofRuns,1.,NofRuns);
4022
4023   Double_t nCMULTot = CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A",
4024                                         seventSelection.Data(),triggerType.Data(),colType.Data())); //Total Nof CMUL7 events
4025
4026   TObjArray* aCluster = striggerCluster.Tokenize("-");
4027   TIter nextCluster(aCluster);
4028   TObjString* striggerClusterS;
4029
4030   TList* lRun2Reject = new TList();
4031   lRun2Reject->SetOwner(kTRUE);
4032
4033   nextRun.Reset();
4034   Int_t i(0);//Run label index
4035   std::cout << std::endl;
4036   std::cout << std::endl;
4037   std::cout << std::endl;
4038   std::cout << "Computing FNorm" << std::endl;
4039   while ( ( s = static_cast<TObjString*>(nextRun())) ) //Run loop
4040   {
4041     Double_t nCMSL(0.),nCMSLandOMUL(0.);
4042     nextCluster.Reset();
4043     while ( (striggerClusterS = static_cast<TObjString*>(nextCluster())) && nCMSL == 0. )
4044     {
4045       nCMSL = CC()->GetSum(Form("/event:%s/trigger:CMSL%s-%s-NOPF-%s/centrality:V0A/run:%s",
4046                                 seventSelection.Data(),triggerType.Data(),colType.Data(),striggerClusterS->GetName(),s->GetName()));
4047
4048       nCMSLandOMUL = CC()->GetSum(Form("/event:%s/trigger:CMSL%s-%s-NOPF-%s&0MUL/centrality:V0A/run:%s",
4049                                        seventSelection.Data(),triggerType.Data(),colType.Data(),striggerClusterS->GetName(),s->GetName()));
4050     }
4051
4052     Double_t nCINT = CC()->GetSum(Form("/event:%s/trigger:CINT%s-%s-NOPF-ALLNOTRD/centrality:V0A/run:%s",
4053                                        seventSelection.Data(),triggerType.Data(),colType.Data(),s->GetName()));
4054
4055     Double_t nCINTandOMSL = CC()->GetSum(Form("/event:%s/trigger:CINT%s-%s-NOPF-ALLNOTRD&0MSL/centrality:V0A/run:%s",
4056                                               seventSelection.Data(),triggerType.Data(),colType.Data(),s->GetName()));
4057
4058     Double_t nCMUL = CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A/run:%s",
4059                                        seventSelection.Data(),triggerType.Data(),colType.Data(),s->GetName()));
4060
4061     Double_t FNorm(0.),FNormError(0.);
4062     Double_t pUfactor = 1.;
4063     if ( nCMSLandOMUL != 0. && nCINTandOMSL !=0. && nCMSL != 0. && nCINT !=0. )
4064     {
4065       if (corrPU)
4066       {
4067         TParameter<Double_t>* p = static_cast<TParameter<Double_t>*>(pUCorr->FindObject(s->GetName()));
4068         if ( p ) pUfactor = p->GetVal();
4069         else
4070         {
4071           AliError(Form("Run %s not found in pile-up correction list",s->GetName()));
4072         }
4073       }
4074       FNorm = (nCMSL*nCINT)*pUfactor/(nCMSLandOMUL*nCINTandOMSL);
4075       FNormError = ErrorPropagationAxBoverCxD(nCMSL,nCINT,nCMSLandOMUL,nCINTandOMSL)*pUfactor;
4076     }
4077     else
4078     {
4079       if ( nCINT == 0 ) std::cout << " Warning: Bad run " << s->GetName() << " has no MB trigger in this bin. Remove from analysis" << std::endl;
4080
4081       lRun2Reject->Add(new TObjString(s->GetName()));
4082       if ( printout ) std::cout << "Run " << s->GetName() << " not used for FNorm cause lack of stats" << std::endl;
4083       continue;
4084     }
4085
4086     FNormTot += FNorm*nCMUL; // This is the sum of equivalent Nof MB per CMUL run by run. NOTE: This sum is NOT always the total equivalent Nof MB per CMUL because in pp 2012 if just one cluster is used at a time this sum is not the sum for all runs
4087     FNormTotError += TMath::Power(nCMUL*FNormError,2.) + TMath::Power(FNorm*TMath::Sqrt(nCMUL),2.);
4088
4089     if ( printout ) std::cout << "Run " << s->GetName() << " FNorm = " << FNorm << " +- " << FNormError << " ; PUFactor =" << pUfactor << " ; " << "Nof CMUL = " << nCMUL << std::endl;
4090
4091     h->GetXaxis()->SetBinLabel(++i,s->GetName());
4092     h->SetBinContent(i,FNorm);
4093     h->SetBinError(i,FNormError);
4094
4095   }
4096
4097   TIter nextRejectRun(lRun2Reject);
4098   TObjString* run2Rej;
4099   Double_t nCMULTotRej(0.);
4100   while ( (run2Rej = static_cast<TObjString*>(nextRejectRun())) )
4101   {
4102     nCMULTotRej += CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A/run:%s",
4103                                      seventSelection.Data(),triggerType.Data(),colType.Data(),
4104                                      run2Rej->GetName())); //Sum of CMUL7 events from rejected runs
4105   }
4106
4107   nCMULTot = nCMULTot - nCMULTotRej;
4108
4109   TH1* o = fMergeableCollection->Histo(id.Data(),h->GetName());
4110
4111   if (o)
4112   {
4113     AliWarning(Form("Replacing %s/%s",id.Data(),h->GetName()));
4114     fMergeableCollection->Remove(Form("%s/%s",id.Data(),h->GetName()));
4115   }
4116
4117   Bool_t adoptOK = fMergeableCollection->Adopt(id.Data(),h);
4118
4119   if ( adoptOK ) std::cout << "+++FNorm histo " << h->GetName() << " adopted" << std::endl;
4120   else AliError(Form("Could not adopt FNorm histo %s",h->GetName()));
4121
4122   //___
4123
4124   FNormTotError =  TMath::Sqrt(TMath::Power(TMath::Sqrt(FNormTotError)/nCMULTot,2.) + TMath::Power(FNormTot*TMath::Sqrt(nCMULTot)/TMath::Power(nCMULTot,2.),2.));
4125
4126   FNormTot = FNormTot/nCMULTot; // nCMULTot is here nCMULTot - nCMULTotRej
4127
4128   std::cout << "FNorm = " << FNormTot << " +- " << FNormTotError << std::endl;
4129
4130   TH1* hFNormTot = new TH1F("hFNormInt","Global Normalization factor",1,0.,1.);
4131   hFNormTot->SetBinContent(1,FNormTot);
4132   hFNormTot->SetBinError(1,FNormTotError);
4133
4134   o = fMergeableCollection->Histo(id.Data(),hFNormTot->GetName());
4135
4136   if (o)
4137   {
4138     AliWarning(Form("Replacing %s/%s",id.Data(),hFNormTot->GetName()));
4139     fMergeableCollection->Remove(Form("%s/%s",id.Data(),hFNormTot->GetName()));
4140   }
4141
4142   adoptOK = fMergeableCollection->Adopt(id.Data(),hFNormTot);
4143
4144   if ( adoptOK ) std::cout << "+++FNorm histo " << hFNormTot->GetName() << " adopted" << std::endl;
4145   else AliError(Form("Could not adopt FNorm histo %s",hFNormTot->GetName()));
4146
4147   //___
4148   TH1* hNEqMB = new TH1F("hNEqMB","Equivalent number of MB events per CMUL",1,0.,1.);
4149
4150   nCMULTot = CC()->GetSum(Form("/event:%s/trigger:CMUL%s-%s-NOPF-MUON/centrality:V0A",
4151                                seventSelection.Data(),triggerType.Data(),colType.Data())); //Total Nof CMUL7 events
4152
4153   Double_t nofEqMB = FNormTot*nCMULTot;
4154   Double_t nofEqMBError = TMath::Sqrt( TMath::Power(FNormTotError*nCMULTot,2.) + TMath::Power(FNormTot*TMath::Sqrt(nCMULTot),2.) );
4155
4156   std::cout << "EqMB = " << nofEqMB << " +- " << TMath::Sqrt(nofEqMBError) << std::endl;
4157
4158   hNEqMB->SetBinContent(1,nofEqMB);
4159   hNEqMB->SetBinError(1,nofEqMBError);
4160
4161   o = fMergeableCollection->Histo(id.Data(),hNEqMB->GetName());
4162
4163   if (o)
4164   {
4165     AliWarning(Form("Replacing %s/%s",id.Data(),hNEqMB->GetName()));
4166     fMergeableCollection->Remove(Form("%s/%s",id.Data(),hNEqMB->GetName()));
4167   }
4168
4169   adoptOK = fMergeableCollection->Adopt(id.Data(),hNEqMB);
4170
4171   if ( adoptOK ) std::cout << "+++FNorm histo " << hNEqMB->GetName() << " adopted" << std::endl;
4172   else AliError(Form("Could not adopt FNorm histo %s",hNEqMB->GetName()));
4173
4174   //___
4175
4176   delete runs;
4177   delete lRun2Reject;
4178   delete aCluster;
4179
4180   return;
4181
4182 }
4183
4184 //_____________________________________________________________________________
4185 void AliAnalysisMuMu::PlotYiedWSyst(const char* triggerCluster)
4186 {
4187   TString striggerCluster(triggerCluster);
4188   if ( striggerCluster.Contains("MUON") && !striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON";
4189   else if ( striggerCluster.Contains("ALLNOTRD") && !striggerCluster.Contains("MUON") ) striggerCluster = "ALLNOTRD";
4190   else if ( striggerCluster.Contains("MUON") && striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON-ALLNOTRD";
4191   else
4192   {
4193     AliError("Unknown trigger cluster");
4194     return;
4195   }
4196
4197   TString path(Form("%s/%s/%s/%s",
4198                     First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kFALSE)).Data(),
4199                     First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data(),
4200                     First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,kFALSE)).Data(),
4201                     First(Config()->GetList(AliAnalysisMuMuConfig::kPairSelectionList,kFALSE)).Data()));
4202
4203   TH1* hY = OC()->Histo(Form("/RESULTS-%s/%s",striggerCluster.Data(),path.Data()),"hJPsiYieldVSdNchdEtaRelative");
4204   if ( !hY )
4205   {
4206     AliError("No yield found");
4207     return;
4208   }
4209
4210   TString id(Form("/TESTSYST/%s",path.Data()));
4211
4212   TH1* hYSyst = static_cast<TH1*>(hY->Clone("RelYieldSyst"));
4213   if ( !hYSyst )
4214   {
4215     AliError("No systematic found");
4216     return;
4217   }
4218
4219   TH1* hS = OC()->Histo(id.Data(),"yield_Systematics");
4220
4221   for ( Int_t i = 1 ; i <= hY->GetNbinsX() ; i++ )
4222   {
4223     hYSyst->SetBinError(i,hS->GetBinContent(i)*hY->GetBinContent(i)/100.);
4224   }
4225
4226   hY->Draw();
4227   hYSyst->Draw("same");
4228 }
4229
4230
4231 //_____________________________________________________________________________
4232 void AliAnalysisMuMu::ComputeJpsiYield(AliMergeableCollection* oc, Bool_t relative, const char* fNormType, const char* triggerCluster,
4233                                        const char* whatever, const char* sResName, AliMergeableCollection* ocMBTrigger, Double_t mNTrCorrection)
4234 {
4235   // This method is suppossed to be used from the file with the counters, oc is the AliMergeableCollection of the file with the histograms (if separated, which is better since we do not need the minv,mean pt... analysis in CINT&0MUL... triggers)
4236   // ocMBTrigger is the mergeableCollection with the MB trigger dNchdEta plot (migth be the same as oc, in which case we set ocMBTrigger=0x0)
4237   //FIXME::Make it general
4238
4239   TString sfNormType(fNormType);
4240   TString swhat("");
4241   TString sres("");
4242   TString swhatever(whatever);
4243   if ( swhatever.Contains("DNCHDETA"))
4244   {
4245     swhat = "dNchdEta";
4246     if ( strlen(sResName) > 0 ) sres = sResName; //"PSIPSIPRIMECB2VWGINDEPTAILS";
4247   }
4248   else if ( swhatever.Contains("NTRCORR") )
4249   {
4250     swhat = "Nch";
4251     if ( strlen(sResName) > 0 ) sres = sResName; //"PSIPSIPRIMECB2VWG_2.0_5.0";
4252   }
4253
4254   if ( IsSimulation() )
4255   {
4256     AliError("Cannot compute J/Psi yield: Is a simulation file");
4257     return;
4258   }
4259
4260   TString striggerCluster(triggerCluster);
4261   if ( striggerCluster.Contains("MUON") && !striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON";
4262   else if ( striggerCluster.Contains("ALLNOTRD") && !striggerCluster.Contains("MUON") ) striggerCluster = "ALLNOTRD";
4263   else if ( striggerCluster.Contains("MUON") && striggerCluster.Contains("ALLNOTRD") ) striggerCluster = "MUON-ALLNOTRD";
4264   else
4265   {
4266     AliError("Unknown trigger cluster");
4267     return;
4268   }
4269
4270   TString path(Form("%s/%s/%s/%s",
4271                     First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kFALSE)).Data(),
4272                     First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data(),
4273                     First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,kFALSE)).Data(),
4274                     First(Config()->GetList(AliAnalysisMuMuConfig::kPairSelectionList,kFALSE)).Data()));
4275
4276   AliMergeableCollection* mc;
4277   if ( !oc ) mc = OC();
4278   else mc = oc;
4279
4280   Double_t bR = 0.0593; // BR(JPsi->mu+mu-)
4281   Double_t bRerror = 0.0006 ;
4282
4283   //_________Integrated yield
4284   AliAnalysisMuMuSpectra* sInt = static_cast<AliAnalysisMuMuSpectra*>(mc->GetObject(Form("/%s/%s",path.Data(),"PSI-INTEGRATED-AccEffCorr")));
4285   if ( !sInt )
4286   {
4287     AliError(Form("No spectra %s found in %s","PSI-INTEGRATED-AccEffCorr",path.Data()));
4288     return;
4289   }
4290
4291   AliAnalysisMuMuBinning* b = new AliAnalysisMuMuBinning;
4292   b->AddBin("psi","INTEGRATED");
4293
4294   AliAnalysisMuMuBinning::Range* bin = static_cast<AliAnalysisMuMuBinning::Range*>(b->CreateBinObjArray()->At(0));
4295
4296   AliAnalysisMuMuResult* result = sInt->GetResultForBin(*bin);
4297   if ( !result )
4298   {
4299     AliError(Form("No result for bin %s found in %s",bin->AsString().Data(),"PSI-INTEGRATED-AccEffCorr"));
4300     return;
4301   }
4302
4303 //  if ( strlen(sResName) > 0/*sResName.Sizeof() > 0*/ )
4304 //  {
4305 //    result = result->SubResult(sres.Data());//INDEPTAILS
4306 //    if ( !result )
4307 //    {
4308 //      AliError(Form("No subresult %s found in %s",sres.Data(),path.Data()));
4309 //      return;
4310 //    }
4311 //  }
4312
4313   Double_t NofJPsiTot = result->GetValue("NofJPsi");
4314   Double_t NofJPsiTotError = result->GetErrorStat("NofJPsi");
4315
4316   TH1* hMBTot = OC()->Histo(Form("/FNORM-%s/PSALL/V0A/hNEqMB",striggerCluster.Data()));//HASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00
4317   if ( !hMBTot )
4318   {
4319     AliError(Form("No eq Nof MB events found in %s",Form("/FNORM-%s/PSALLHASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00/V0A/hNEqMB",striggerCluster.Data())));
4320     return;
4321   }
4322
4323   Double_t nEqMBTot = hMBTot->GetBinContent(1);
4324   Double_t nEqMBTotError = hMBTot->GetBinError(1);
4325
4326   Double_t yieldInt = NofJPsiTot/(nEqMBTot*bR);
4327   Double_t yieldIntError = TMath::Sqrt(TMath::Power(NofJPsiTotError/(nEqMBTot*bR),2.) +
4328                                        TMath::Power(nEqMBTotError*NofJPsiTot*bR/TMath::Power(nEqMBTot*bR,2.),2.) +
4329                                        TMath::Power(NofJPsiTot*nEqMBTot*bRerror/TMath::Power(nEqMBTot*bR,2.),2.));
4330
4331   std::cout << "Integrated yield = " << yieldInt << " +- " << yieldIntError << std::endl;
4332
4333   TH1* hYint = new TH1F("hJPsiYieldInt","Integrated J/#psi yield",1,0.,1.);
4334   hYint->SetBinContent(1,yieldInt);
4335   hYint->SetBinError(1,yieldIntError);
4336
4337   TH1* o = mc->Histo(Form("/RESULTS-%s/%s",striggerCluster.Data(),path.Data()),hYint->GetName());
4338
4339   if (o)
4340   {
4341     AliWarning(Form("Replacing /RESULTS-%s/%s/%s",striggerCluster.Data(),path.Data(),hYint->GetName()));
4342     mc->Remove(Form("/RESULTS-%s/%s/%s",striggerCluster.Data(),path.Data(),hYint->GetName()));
4343   }
4344
4345   Bool_t adoptOK = mc->Adopt(Form("/RESULTS-%s/%s",striggerCluster.Data(),path.Data()),hYint);
4346
4347   if ( adoptOK ) std::cout << "+++Yield histo " << hYint->GetName() << " adopted" << std::endl;
4348   else AliError(Form("Could not adopt Yield histo %s",hYint->GetName()));
4349
4350   delete b;
4351
4352   //_____Differential yield
4353
4354   AliAnalysisMuMuSpectra* s = static_cast<AliAnalysisMuMuSpectra*>(mc->GetObject(Form("/%s/%s",path.Data(),whatever)));
4355   if ( !s )
4356   {
4357     AliError(Form("No spectra %s found in %s",whatever,path.Data()));
4358     return;
4359   }
4360
4361   std::cout << "Number of J/Psi:" << std::endl;
4362   TH1* hry = s->Plot("NofJPsi",sres.Data(),kFALSE);//INDEPTAILS //Number of Jpsi
4363
4364   std::cout << "" << std::endl;
4365
4366 //  std::cout << "Equivalent number of MB events:" << std::endl;
4367   TH1* hMB(0x0);
4368   if ( sfNormType.Contains("offline") )
4369   {
4370     hMB = OC()->Histo(Form("/FNORM-%s/PSALL/V0A/hNofEqMBVSdNchdEta",striggerCluster.Data()));//HASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00
4371     if ( !hMB )
4372     {
4373       AliError("Histo hNofEqMBVSdNchdEta not found");
4374       return;
4375     }
4376   }
4377   else if ( sfNormType.Contains("global") )
4378   {
4379     hMB = OC()->Histo(Form("/FNORM-%s/PSALL/V0A/hNofEqMBVSdNchdEtaFromGlobal",striggerCluster.Data()));//HASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00
4380     if ( !hMB )
4381     {
4382       AliError("Histo hNofEqMBVSdNchdEtaFromGlobal not found");
4383       return;
4384     }
4385   }
4386   else if ( sfNormType.Contains("mean") )
4387   {
4388     hMB = OC()->Histo(Form("/FNORM-%s/PSALL/V0A/hNofEqMBVSdNchdEtaFromMean",striggerCluster.Data()));//HASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00
4389     if ( !hMB )
4390     {
4391       AliError("Histo hNofEqMBVSdNchdEtaFromMean not found");
4392       return;
4393     }
4394   }
4395   else
4396   {
4397     AliError("Dont know what Fnorm use");
4398     return;
4399   }
4400
4401   TH1* hy;
4402   if ( relative )
4403   {
4404     TString path2(Form("/%s/%s/%s",
4405                       First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kFALSE)).Data(),
4406                       First(Config()->GetList(AliAnalysisMuMuConfig::kMinbiasTriggerList,kFALSE)).Data(),
4407                       First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,kFALSE)).Data()));
4408
4409     TH1* hdNch;
4410     if ( ocMBTrigger ) hdNch = ocMBTrigger->Histo(path2.Data(),swhat.Data());//dNchdEta
4411     else hdNch = mc->Histo(path2.Data(),swhat.Data());//dNchdEta
4412
4413     const TArrayD* binArray = hry->GetXaxis()->GetXbins();
4414     Int_t size = binArray->GetSize();
4415     Double_t* axis = new Double_t[size];
4416     for ( Int_t k = 0 ; k < size ; k++ )
4417     {
4418       axis[k] = binArray->At(k)/(hdNch->GetMean()*(1 - mNTrCorrection));
4419     }
4420
4421     hy = new TH1D("hJPsiYieldVSdNchdEtaRelative","Relative J/#psi yield vs dN_{ch}/d#eta;dN_{ch}/d#eta/<dN_{ch}/d#eta>;Y^{J/#psi}/Y^{J/#psi}_{int}",size-1,axis);
4422     delete axis;
4423   }
4424   else
4425   {
4426     hy = static_cast<TH1D*>(hry->Clone("hJPsiYieldVSdNchdEta"));
4427     hy->SetTitle("J/#psi yield vs dN_{ch}/d#eta");
4428     hy->GetXaxis()->SetTitle("dN_{ch}/d#eta");
4429     hy->GetYaxis()->SetTitle("Y^{J/#psi}");
4430   }
4431                                     // AccxEff(from rel diff or paper)  // Signal extraction
4432 //  Double_t systNofJpsiBin[9] = {TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.01,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.008,2.) ),TMath::Sqrt( TMath::Power(0.022,2.) + TMath::Power(0.007,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.008,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.007,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.009,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.008,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.016 ,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.033,2.) )}; //FIXME: find a way to give this as input
4433 //  Double_t systFNorm[9] = {0.003,0.001,0.002,0.003,0.002,0.004,0.011,0.012,0.071};
4434 //  Double_t systPU[9] = {0.00,0.01,0.012,0.014,0.014,0.019,0.020,0.021,0.040}; //_______pPb
4435        // AccxEff(from paper)
4436 //  Double_t systNofJpsiTot = 0.015;
4437
4438 //  Double_t systNofJpsiBin[9] = {TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.007,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.006,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.005,2.) ),TMath::Sqrt( TMath::Power(0.028,2.) + TMath::Power(0.006,2.) ),TMath::Sqrt( TMath::Power(0.016,2.) + TMath::Power(0.004,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.004,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.006,2.) ),TMath::Sqrt( TMath::Power(0.024,2.) + TMath::Power(0.005 ,2.) ),TMath::Sqrt( TMath::Power(0.015,2.) + TMath::Power(0.016,2.) )}; //FIXME: find a way to give this as input
4439 //  Double_t systFNorm[9] = {0.005,0.004,0.004,0.004,0.003,0.002,0.002,0.04,0.04};
4440 //  Double_t systPU[9] = {0.00,0.007,0.015,0.011,0.014,0.018,0.014,0.011,0.020}; //______Pbp
4441 //  Double_t systNofJpsiTot = 0.015;
4442
4443 //  Double_t systNofJpsiBin[9] = {TMath::Sqrt( TMath::Power(0.034,2.) + TMath::Power(0.005,2.) ),TMath::Sqrt( TMath::Power(0.017,2.) + TMath::Power(0.005,2.) ),TMath::Sqrt( TMath::Power(0.017,2.) + TMath::Power(0.004,2.) ),TMath::Sqrt( TMath::Power(0.017,2.) + TMath::Power(0.005,2.) ),TMath::Sqrt( TMath::Power(0.042,2.) + TMath::Power(0.002,2.) ),TMath::Sqrt( TMath::Power(0.063,2.) + TMath::Power(0.014,2.) ),TMath::Sqrt( TMath::Power(0.094,2.) + TMath::Power(0.009,2.) ),TMath::Sqrt( TMath::Power(0.00,2.) + TMath::Power(0.00 ,2.) ),TMath::Sqrt( TMath::Power(0.00,2.) + TMath::Power(0.00,2.) )}; //FIXME: find a way to give this as input
4444 //  Double_t systFNorm[9] = {0.004,0.019,0.002,0.012,0.048,0.063,0.082,0.000,0.000};
4445 //  Double_t systPU[9] = {0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00}; //______pp |eta|<0.5
4446 //  Double_t systNofJpsiTot = 0.017;
4447
4448   Double_t systNofJpsiBin[9] = {TMath::Sqrt( TMath::Power(0.037,2.) + TMath::Power(0.002,2.) ),TMath::Sqrt( TMath::Power(0.021,2.) + TMath::Power(0.002,2.) ),TMath::Sqrt( TMath::Power(0.022,2.) + TMath::Power(0.002,2.) ),TMath::Sqrt( TMath::Power(0.017,2.) + TMath::Power(0.002,2.) ),TMath::Sqrt( TMath::Power(0.019,2.) + TMath::Power(0.001,2.) ),TMath::Sqrt( TMath::Power(0.036,2.) + TMath::Power(0.002,2.) ),TMath::Sqrt( TMath::Power(0.042,2.) + TMath::Power(0.001,2.) ),TMath::Sqrt( TMath::Power(0.039,2.) + TMath::Power(0.012 ,2.) ),TMath::Sqrt( TMath::Power(0.000,2.) + TMath::Power(0.000,2.) )}; //FIXME: find a way to give this as input
4449   Double_t systFNorm[9] = {0.026,0.002,0.015,0.019,0.012,0.030,0.015,0.119,0.000};
4450   Double_t systPU[9] = {0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00}; //______pp |eta|<1
4451   Double_t systNofJpsiTot = 0.017;
4452
4453   for ( Int_t i = 1 ; i <= hy->GetNbinsX() ; i++ )
4454   {
4455     Double_t yield = hry->GetBinContent(i)/(hMB->GetBinContent(i)*bR);
4456     Double_t yieldError = TMath::Sqrt(TMath::Power(hry->GetBinError(i)/(hMB->GetBinContent(i)*bR),2.) +
4457                                       TMath::Power(hMB->GetBinError(i)*hry->GetBinContent(i)*bR/TMath::Power(hMB->GetBinContent(i)*bR,2.),2.) +
4458                                       TMath::Power(hry->GetBinContent(i)*hMB->GetBinContent(i)*bRerror/TMath::Power(hMB->GetBinContent(i)*bR,2.),2.));
4459
4460 //    std::cout << "Differential yield bin " << i << " = " << yield << " +- " << yieldError << std::endl;
4461
4462     if ( relative )
4463     {
4464       yieldError = TMath::Sqrt(TMath::Power(yieldError/yieldInt,2.) + TMath::Power((yield*yieldIntError)/TMath::Power(yieldInt,2.),2.));
4465       yield /= yieldInt;
4466
4467 //      std::cout << "relative yield bin " << i << " = " << yield << " +- " << yieldError << std::endl;
4468       Double_t sNJpsiBin = hry->GetBinContent(i)*systNofJpsiBin[i-1];
4469       Double_t sNJpsiTot = NofJPsiTot*systNofJpsiTot;
4470       Double_t sMBBin = hMB->GetBinContent(i)*systFNorm[i-1];
4471       Double_t sMBTot = nEqMBTot*0.01;
4472
4473       Double_t syst = TMath::Sqrt( TMath::Power((sNJpsiBin/NofJPsiTot)*(nEqMBTot/hMB->GetBinContent(i)),2.) + TMath::Power((hry->GetBinContent(i)*sNJpsiTot/TMath::Power(NofJPsiTot,2.))*(nEqMBTot/hMB->GetBinContent(i)),2.) + TMath::Power((hry->GetBinContent(i)/NofJPsiTot)*(sMBTot/hMB->GetBinContent(i)),2.) + TMath::Power((hry->GetBinContent(i)/NofJPsiTot)*(sMBBin*nEqMBTot/TMath::Power(hMB->GetBinContent(i),2.)),2.) );
4474
4475       std::cout << "sys" << syst/yield << " w/pu = " << TMath::Sqrt( TMath::Power(syst/yield,2.) + TMath::Power(systPU[i-1],2.)) << std::endl;
4476       std::cout << yield << " +- " << yieldError << std::endl;
4477     }
4478
4479     hy->SetBinContent(i,yield);
4480     hy->SetBinError(i,yieldError);
4481   }
4482
4483   o = mc->Histo(Form("/RESULTS-%s/%s",striggerCluster.Data(),path.Data()),hy->GetName());
4484
4485   if (o)
4486   {
4487     AliWarning(Form("Replacing %s/%s","/RESULTS-%s/PSALLHASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00/V0A",hy->GetName()));
4488     mc->Remove(Form("/RESULTS-%s/%s/%s",striggerCluster.Data(),path.Data(),hy->GetName()));
4489   }
4490
4491   adoptOK = mc->Adopt(Form("/RESULTS-%s/%s",striggerCluster.Data(),path.Data()),hy);
4492
4493   if ( adoptOK ) std::cout << "+++Yield histo " << hy->GetName() << " adopted" << std::endl;
4494   else AliError(Form("Could not adopt Yield histo %s",hy->GetName()));
4495
4496
4497
4498   delete hry;
4499
4500
4501   return;
4502 }
4503
4504 //_____________________________________________________________________________
4505 void AliAnalysisMuMu::ComputeJpsiMPt(Bool_t relative, const char* whatever, const char* sResName, AliMergeableCollection* ocMBTrigger, Double_t mNTrCorrection)
4506 {
4507   // ocMBTrigger is the mergeableCollection with the MB trigger dNchdEta plot (migth be the same as oc, in which case we set ocMBTrigger=0x0)
4508   //FIXME::Make it general
4509
4510
4511   TString swhat("");
4512   TString sres("");
4513   TString swhatever(whatever);
4514 //  if ( swhatever.Contains("DNCHDETA"))
4515 //  {
4516 //    swhat = "dNchdEta";
4517 //    sres = "MPT2CB2VWGPOL2INDEPTAILS";
4518 //  }
4519 //  else
4520     if ( swhatever.Contains("NTRCORR") )
4521   {
4522     swhat = "Nch";
4523     if ( strlen(sResName) > 0 ) sres = sResName; //sres = "MPTPSIPSIPRIMECB2VWG_BKGMPTPOL2";
4524   }
4525
4526   if ( IsSimulation() )
4527   {
4528     AliError("Cannot compute J/Psi yield: Is a simulation file");
4529     return;
4530   }
4531
4532   TString path(Form("%s/%s/%s/%s",
4533                     First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kFALSE)).Data(),
4534                     First(Config()->GetList(AliAnalysisMuMuConfig::kDimuonTriggerList,kFALSE)).Data(),
4535                     First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,kFALSE)).Data(),
4536                     First(Config()->GetList(AliAnalysisMuMuConfig::kPairSelectionList,kFALSE)).Data()));
4537
4538   //_________Integrated mean pt
4539   AliAnalysisMuMuSpectra* sInt = static_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(Form("/%s/%s",path.Data(),"PSI-INTEGRATED-AccEffCorr-MeanPtVsMinvUS")));
4540   if ( !sInt )
4541   {
4542     AliError(Form("No spectra %s found in %s","PSI-INTEGRATED-AccEffCorr-MeanPtVsMinvUS",path.Data()));
4543     return;
4544   }
4545
4546   AliAnalysisMuMuBinning* b = new AliAnalysisMuMuBinning;
4547   b->AddBin("psi","INTEGRATED");
4548
4549   AliAnalysisMuMuBinning::Range* bin = static_cast<AliAnalysisMuMuBinning::Range*>(b->CreateBinObjArray()->At(0));
4550
4551   AliAnalysisMuMuResult* result = sInt->GetResultForBin(*bin);
4552   if ( !result )
4553   {
4554     AliError(Form("No result for bin %s found in spectra %s",bin->AsString().Data(),sInt->GetName()));
4555     return;
4556   }
4557
4558 //  if ( sres.Sizeof() > 0 )
4559 //  {
4560 //    result = result->SubResult(sres.Data());
4561 //    if ( !result )
4562 //    {
4563 //      AliError(Form("No subresult %s found in result",result->GetName()));
4564 //      return;
4565 //    }
4566 //  AliAnalysisMuMuResult* subresult = result->SubResult(sres.Data());//"MPT2CB2VWGPOL2INDEPTAILS"
4567 //  if ( !subresult )
4568 //  {
4569 //    AliError(Form("No subresult MPT2CB2VWGPOL2 found in result %s",result->GetName()));
4570 //    return;
4571 //  }
4572
4573 //  }
4574 //  Double_t JPsiMPtTot = subresult->GetValue("MeanPtJPsi");
4575 //  Double_t JPsiMPtTotError = subresult->GetErrorStat("MeanPtJPsi");
4576
4577   Double_t JPsiMPtTot = result->GetValue("MeanPtJPsi");
4578   Double_t JPsiMPtTotError = result->GetErrorStat("MeanPtJPsi");
4579
4580   TH1* hMPtint = new TH1F("hJPsiMPtInt","Integrated J/#psi mean p_{T}",1,0.,1.);
4581   hMPtint->SetBinContent(1,JPsiMPtTot);
4582   hMPtint->SetBinError(1,JPsiMPtTotError);
4583
4584   TH1* o = OC()->Histo(Form("/RESULTS/%s",path.Data()),hMPtint->GetName());
4585
4586   if (o)
4587   {
4588     AliWarning(Form("Replacing /RESULTS/%s/%s",path.Data(),hMPtint->GetName()));
4589     OC()->Remove(Form("/RESULTS/%s/%s",path.Data(),hMPtint->GetName()));
4590   }
4591
4592   Bool_t adoptOK = OC()->Adopt(Form("/RESULTS/%s",path.Data()),hMPtint);
4593
4594   if ( adoptOK ) std::cout << "+++Mean Pt histo " << hMPtint->GetName() << " adopted" << std::endl;
4595   else AliError(Form("Could not adopt Mean Pt histo %s",hMPtint->GetName()));
4596
4597   delete b;
4598
4599    //_____Differential mean pt
4600
4601   AliAnalysisMuMuSpectra* s = static_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(Form("/%s/%s",path.Data(),whatever)));
4602   if ( !s )
4603   {
4604     AliError(Form("No spectra %s found in %s",whatever,path.Data()));
4605     return;
4606   }
4607
4608   std::cout << "Mean pt of J/Psi:" << std::endl;
4609   TH1* hrmPt = s->Plot("MeanPtJPsi",sres.Data(),kFALSE); //MPT2CB2VWGPOL2INDEPTAILS//mean pt of Jpsi
4610   std::cout << "" << std::endl;
4611
4612   Double_t ptInt,ptIntError;
4613   TH1* hmPt;
4614   if ( relative )
4615   {
4616     TString path2(Form("/%s/%s/%s",
4617                        First(Config()->GetList(AliAnalysisMuMuConfig::kEventSelectionList,kFALSE)).Data(),
4618                        First(Config()->GetList(AliAnalysisMuMuConfig::kMinbiasTriggerList,kFALSE)).Data(),
4619                        First(Config()->GetList(AliAnalysisMuMuConfig::kCentralitySelectionList,kFALSE)).Data()));
4620
4621     TH1* hdNch;
4622     if ( ocMBTrigger ) hdNch = ocMBTrigger->Histo(path2.Data(),swhat.Data());
4623     else hdNch = OC()->Histo(path2.Data(),swhat.Data());
4624
4625     const TArrayD* binArray = hrmPt->GetXaxis()->GetXbins();
4626     Int_t size = binArray->GetSize();
4627     Double_t* axis = new Double_t[size];
4628     for ( Int_t k = 0 ; k < size ; k++ )
4629     {
4630       axis[k] = binArray->At(k)/(hdNch->GetMean()*(1 - mNTrCorrection));
4631     }
4632
4633     hmPt = new TH1D("hJPsiMeanPtVSdNchdEtaRelative","Relative J/#psi mean p_{T} vs dN_{ch}/d#eta/<dN_{ch}/d#eta>;dN_{ch}/d#eta/<dN_{ch}/d#eta>;<p_{T}^{J/#psi}>/<p_{T}^{J/#psi}_{int}>",size-1,axis);
4634     delete axis;
4635
4636     ptInt = result->GetValue("MeanPtJPsi",sres.Data());
4637     ptIntError = result->GetErrorStat("MeanPtJPsi",sres.Data());
4638
4639 //    delete b;
4640   }
4641   else
4642   {
4643     hmPt = static_cast<TH1D*>(hrmPt->Clone("hJPsiMeanPtVSdNchdEta"));
4644     hmPt->SetTitle("J/#psi mean p_{T} vs dN_{ch}/d#eta");
4645     hmPt->GetXaxis()->SetTitle("dN_{ch}/d#eta");
4646     hmPt->GetYaxis()->SetTitle("<p_{T}^{J/#psi}>");
4647   }
4648
4649   Double_t systMptInt[9] = {0.014,0.014,0.014,0.014,0.014,0.014,0.014,0.014,0.014}; //FIXME: find a way to give this as input
4650   Double_t systMptBin[9] = {0.014,0.014,0.014,0.014,0.014,0.014,0.014,0.014,0.014};
4651
4652   Double_t systMptRel[9] = {0.002,0.001,0.001,0.002,0.002,0.002,0.002,0.004,0.004}; //signal extraction pPb
4653
4654 //  Double_t systMptRel[9] = {0.002,0.001,0.012,0.001,0.002,0.002,0.001,0.004,0.003}; //signal extraction Pbp
4655
4656 //  Double_t systMptRel[9] = {0.001,0.002,0.001,0.002,0.002,0.003,0.005,0.000,0.000}; //signal extraction pp|eta|<05
4657
4658 //  Double_t systMptRel[9] = {0.002,0.002,0.002,0.002,0.001,0.002,0.001,0.012,0.000}; //signal extraction pp|eta|<1
4659
4660   for ( Int_t i = 1 ; i <= hrmPt->GetNbinsX() ; i++ )
4661   {
4662     Double_t pt = hrmPt->GetBinContent(i);
4663     Double_t ptError = hrmPt->GetBinError(i);
4664
4665     if ( relative )
4666     {
4667       ptError = TMath::Sqrt(TMath::Power(ptError/ptInt,2.) + TMath::Power((pt*ptIntError)/TMath::Power(ptInt,2.),2.));
4668
4669       Double_t sMptInt = ptInt*systMptInt[i-1];
4670       Double_t sMptBin = pt*systMptBin[i-1];
4671       Double_t sysMptRel = TMath::Sqrt( TMath::Power(sMptBin/ptInt,2) + TMath::Power(pt*sMptInt/TMath::Power(ptInt,2.),2.) );
4672
4673       pt /= ptInt;
4674
4675       std::cout << TMath::Sqrt( TMath::Power(sysMptRel/pt,2.) +TMath::Power(systMptRel[i-1],2.) ) << std::endl;
4676
4677       std::cout << pt << " +- " << ptError << std::endl;
4678
4679     }
4680
4681     hmPt->SetBinContent(i,pt);
4682     hmPt->SetBinError(i,ptError);
4683   }
4684
4685   o = fMergeableCollection->Histo(Form("/RESULTS/%s",path.Data()),hmPt->GetName());
4686
4687   if (o)
4688   {
4689     AliWarning(Form("Replacing /RESULTS/%s/%s",path.Data(),hmPt->GetName()));
4690     fMergeableCollection->Remove(Form("/RESULTS/%s/%s",path.Data(),hmPt->GetName()));
4691   }
4692
4693   adoptOK = fMergeableCollection->Adopt(Form("/RESULTS/%s",path.Data()),hmPt);
4694
4695   if ( adoptOK ) std::cout << "+++Mean Pt histo " << hmPt->GetName() << " adopted" << std::endl;
4696   else AliError(Form("Could not adopt mean pt histo %s",hmPt->GetName()));
4697
4698
4699
4700   delete hrmPt;
4701
4702
4703   return;
4704
4705
4706 }
4707
4708 //_____________________________________________________________________________
4709 Double_t AliAnalysisMuMu::ErrorPropagationAxBoverCxD(Double_t a,Double_t b,Double_t c,Double_t d)
4710 {
4711   //Just valid for counts
4712   Double_t error2 = TMath::Power(b/(c*d),2.)*a + TMath::Power(a/(c*d),2.)*b + TMath::Power(a*b*d,2.)*(c/TMath::Power(c*d,4.)) + TMath::Power(a*b*c,2.)*(d/TMath::Power(c*d,4.));
4713
4714   return TMath::Sqrt(error2);
4715 }
4716
4717 //_____________________________________________________________________________
4718 TH1* AliAnalysisMuMu::ComputeEquNofMB(const char* what,const char* quantity,const char* flavour,Bool_t printout)
4719 {
4720
4721   AliAnalysisMuMuBinning* binning = BIN()->Project(what,quantity,flavour);
4722   TObjArray* dNchdEtas = binning->CreateBinObjArray();
4723
4724   Double_t* binArray = binning->CreateBinArray();
4725
4726   TIter next(dNchdEtas);
4727   AliAnalysisMuMuBinning::Range* r;
4728
4729   TH1* hFNorm = ComputeDiffFnormFromHistos(what,quantity,flavour,kFALSE);
4730
4731   TH1* hNMB = new TH1F("hNofEqMB","Equivalent number of MB triggers vs dN_{ch}/d#eta;dN_{ch}/d#eta;FNorm",dNchdEtas->GetEntries(),binArray);
4732
4733   Int_t bin(0);
4734   while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(next()) ) )
4735   {
4736
4737     TH1* hCMUL = OC()->Histo(Form("/PSALLHASSPDSPDZQA_RES0.25_ZDIF0.50SPDABSZLT10.00/CMUL7-B-NOPF-MUON/V0A/%s",
4738                                   Form("EventsIn%s",r->AsString().Data())));
4739     if ( !hCMUL )
4740     {
4741       AliError(Form("No event histo in bin %s found for CMUL7-B-NOPF-MUON",r->AsString().Data()));
4742       return 0x0;
4743     }
4744
4745     Double_t NMB = hCMUL->GetBinContent(1)*hFNorm->GetBinContent(++bin);
4746     Double_t NMBError = TMath::Sqrt(TMath::Power(hCMUL->GetBinContent(1)*hFNorm->GetBinError(bin),2.) + TMath::Power(TMath::Sqrt(hCMUL->GetBinContent(1))*hFNorm->GetBinContent(bin),2));
4747
4748     if ( printout ) std::cout << r->AsString().Data() << " : " << NMB << " +- " << NMBError << std::endl;
4749
4750     hNMB->SetBinContent(bin,NMB);
4751     hNMB->SetBinError(bin,NMBError);
4752   }
4753
4754   delete dNchdEtas;
4755   delete[] binArray;
4756
4757   return hNMB;
4758 }
4759
4760
4761 //_____________________________________________________________________________
4762 AliAnalysisMuMuSpectra* AliAnalysisMuMu::CorrectSpectra(const char* type, const char* flavour)
4763 {
4764   /// Correct one spectra
4765
4766   if (!SIM())
4767   {
4768     AliError("Cannot compute corrected yield without associated MC file !");
4769     return 0x0;
4770   }
4771
4772   const char* accEffSubResultName="PSICOUNT:1";
4773
4774   AliAnalysisMuMuSpectra* realSpectra = GetSpectra(type,flavour);
4775   AliAnalysisMuMuSpectra* simSpectra = SIM()->GetSpectra(type,flavour);
4776
4777   if ( !realSpectra )
4778   {
4779     AliError("could not get real spectra");
4780     return 0x0;
4781   }
4782
4783   if ( !simSpectra)
4784   {
4785     AliError("could not get sim spectra");
4786     return 0x0;
4787   }
4788
4789   realSpectra->Correct(*simSpectra,"Jpsi",accEffSubResultName);
4790
4791   Update();
4792
4793   return realSpectra;
4794 }
4795
4796 //_____________________________________________________________________________
4797 AliAnalysisMuMuSpectra* AliAnalysisMuMu::ComputeYield(const char* type, const char* flavour)
4798 {
4799   if (!SIM())
4800   {
4801     AliError("Cannot compute corrected yield without associated MC file !");
4802     return 0x0;
4803   }
4804
4805   const char* accEffSubResultName="PSICOUNT:1";
4806
4807   AliAnalysisMuMuSpectra* realSpectra = GetSpectra(type,flavour);
4808
4809   if ( !realSpectra )
4810   {
4811     AliError("could not get real spectra");
4812     return 0x0;
4813   }
4814
4815   if (!realSpectra->HasValue("CoffNofJpsi"))
4816   {
4817     if (!CorrectSpectra(type,flavour))
4818     {
4819       AliError("Could not get corrected spectra");
4820       return 0x0;
4821     }
4822   }
4823
4824   AliAnalysisMuMuSpectra* simSpectra = SIM()->GetSpectra(type,flavour);
4825
4826   if ( !simSpectra)
4827   {
4828     AliErrorClass("could not get sim spectra");
4829     return 0x0;
4830   }
4831
4832   Double_t nofCMUL7 = CC()->GetSum(Form("trigger:CMUL7-B-NOPF-MUON/event:PSALL"));
4833   Double_t nofCINT7 = CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD/event:PSALL"));
4834   Double_t nofCINT7w0MUL = CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD&0MUL/event:PSALL"));
4835
4836   AliAnalysisMuMuBinning* binning = realSpectra->Binning();
4837   TObjArray* bins = binning->CreateBinObjArray();
4838   TIter nextBin(bins);
4839   AliAnalysisMuMuBinning::Range* bin;
4840   Int_t i(0);
4841   AliAnalysisMuMuJpsiResult* r;
4842
4843   while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
4844   {
4845     r = static_cast<AliAnalysisMuMuJpsiResult*>(realSpectra->BinContentArray()->At(i));
4846
4847     StdoutToAliDebug(1,std::cout << "bin=";r->Print(););
4848
4849     AliAnalysisMuMuJpsiResult* rsim = static_cast<AliAnalysisMuMuJpsiResult*>(simSpectra->BinContentArray()->At(i));
4850
4851     Double_t mbeq = nofCINT7w0MUL / ( nofCINT7 * nofCMUL7);
4852     Double_t mbeqError = mbeq * AliAnalysisMuMuResult::ErrorABC( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
4853                                                                 nofCINT7,TMath::Sqrt(nofCINT7),
4854                                                                 nofCMUL7,TMath::Sqrt(nofCMUL7));
4855
4856     r->Set("Fnorm",nofCINT7/nofCINT7w0MUL,(nofCINT7/nofCINT7w0MUL)*AliAnalysisMuMuResult::ErrorAB( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
4857                                                                                                 nofCINT7,TMath::Sqrt(nofCINT7)));
4858
4859     Double_t yield =  r->GetValue("CorrNofJpsi") * mbeq;
4860
4861     Double_t yieldError = yield * AliAnalysisMuMuResult::ErrorAB( r->GetValue("CorrNofJpsi"), r->GetErrorStat("CorrNofJpsi"),
4862                                                                  mbeq,mbeqError);
4863
4864     r->Set("YJpsi",yield,yieldError);
4865
4866     r->Set("NofInputJpsi",rsim->GetValue("NofInputJpsi",accEffSubResultName),rsim->GetErrorStat("NofInputJpsi",accEffSubResultName));
4867     r->Set("AccEffJpsi",rsim->GetValue("AccEffJpsi",accEffSubResultName),rsim->GetErrorStat("AccEffJpsi",accEffSubResultName));
4868
4869     ++i;
4870   }
4871
4872   delete bins;
4873
4874   Update();
4875
4876   return realSpectra;
4877 }
4878
4879 //_____________________________________________________________________________
4880 AliAnalysisMuMuSpectra* AliAnalysisMuMu::RABy(const char* type, const char* direction)
4881 {
4882   /// Compute the RAB...
4883
4884   if (!SIM()) return 0x0;
4885
4886   Double_t rapidityShift = 0.465;// 0.5*TMath::Log(208.0/82.0);
4887   const Double_t sqrts=5.023;
4888   const Double_t ymax=TMath::Log(sqrts*1000.0/3.096916);
4889   const Double_t tab = 0.093e-6; // nb^-1
4890   const Double_t tabError = 0.0035E-6; // nb^-1
4891   const char* accEffSubResultName="PSICOUNT:1";
4892
4893   TF1 ydist("ydist","[0]*TMath::Exp(-(x*x)/(2.0*0.39*0.39))",0.,0.5);
4894   ydist.SetParameter(0,1.);
4895
4896   //Normalization to the values presented by Zaida and Rosana on January 11th 2013 https://indico.cern.ch/conferenceDisplay.py?confId=224985 slide 22
4897   // Normalization is done in the rapidity range 2.75<y<3.25 where Rosanas values is 230.8+212.1
4898   Double_t y1_norma= 2.75/ymax;
4899   Double_t y2_norma= 3.25/ymax;
4900   Double_t normalization = 0.25*(230.8+212.1)/ydist.Integral(y1_norma, y2_norma);
4901   ydist.SetParameter(0,normalization);
4902 //  AliInfoClass(Form("ymax=%e normalization=%f",ymax,ydist.Integral(y1_norma, y2_norma)));
4903
4904   AliAnalysisMuMuSpectra* realSpectra = static_cast<AliAnalysisMuMuSpectra*>(OC()->GetObject(Form("/PSALL/CMUL7-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-%s",type)));
4905   AliAnalysisMuMuSpectra* simSpectra = static_cast<AliAnalysisMuMuSpectra*>(SIM()->OC()->GetObject(Form("/ALL/CMULLO-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-%s",type)));
4906
4907   if ( !realSpectra )
4908   {
4909     AliErrorClass("could not get real spectra");
4910     return 0x0;
4911   }
4912
4913   if ( !simSpectra)
4914   {
4915     AliErrorClass("could not get sim spectra");
4916     return 0x0;
4917   }
4918
4919   AliAnalysisMuMuSpectra* corrSpectra = static_cast<AliAnalysisMuMuSpectra*>(realSpectra->Clone());
4920   corrSpectra->Correct(*simSpectra,"Jpsi",accEffSubResultName);
4921
4922   Double_t nofCMUL7 = CC()->GetSum(Form("trigger:CMUL7-B-NOPF-MUON/event:PSALL"));
4923   Double_t nofCINT7 = CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD/event:PSALL"));
4924   Double_t nofCINT7w0MUL = CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD&0MUL/event:PSALL"));
4925
4926   AliAnalysisMuMuBinning* binning = realSpectra->Binning();
4927   TObjArray* bins = binning->CreateBinObjArray();
4928   TIter nextBin(bins);
4929   AliAnalysisMuMuBinning::Range* bin;
4930   Int_t i(0);
4931   AliAnalysisMuMuJpsiResult* r;
4932
4933   Int_t n = bins->GetLast();
4934
4935   TObjArray finalBins(n+1);
4936   finalBins.SetOwner(kTRUE);
4937
4938   TObjArray finalResults(n+1);
4939   finalResults.SetOwner(kFALSE);
4940
4941   while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
4942   {
4943     Double_t ylowlab = bin->Xmin();
4944     Double_t yhighlab = bin->Xmax();
4945
4946     Double_t ylowcms, yhighcms;
4947     Double_t ylownorm, yhighnorm;
4948
4949     if ( bin->IsIntegrated() )
4950     {
4951       ylowlab = -4;
4952       yhighlab = -2.5;
4953     }
4954
4955     if ( strcmp(direction,"pPb")==0 )
4956     {
4957       ylowcms = TMath::Abs(yhighlab) -  rapidityShift;
4958       yhighcms = TMath::Abs(ylowlab) - rapidityShift;
4959       ylownorm = ylowcms/ymax;
4960       yhighnorm = yhighcms/ymax;
4961     }
4962     else
4963     {
4964       ylowcms = ylowlab - rapidityShift;
4965       yhighcms = yhighlab - rapidityShift;
4966       ylownorm = -yhighcms/ymax;
4967       yhighnorm = -ylowcms/ymax;
4968     }
4969
4970
4971     Double_t brsigmapp = ydist.Integral(ylownorm,yhighnorm);
4972     Double_t brsigmappError = 0.0; // FIXME
4973
4974     AliInfoClass(Form("y range : LAB %f ; %f CMS %f ; %f -> ynorm : %f ; %f -> BR x sigmapp = %f",
4975                       ylowlab,yhighlab,ylowcms,yhighcms,ylownorm,yhighnorm,brsigmapp));
4976
4977     r = static_cast<AliAnalysisMuMuJpsiResult*>(corrSpectra->BinContentArray()->At(i)->Clone());
4978
4979     AliAnalysisMuMuJpsiResult* rsim = static_cast<AliAnalysisMuMuJpsiResult*>(simSpectra->BinContentArray()->At(i));
4980
4981     Double_t mbeq = nofCINT7w0MUL / ( nofCINT7 * nofCMUL7);
4982     Double_t mbeqError = mbeq * AliAnalysisMuMuResult::ErrorABC( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
4983                                          nofCINT7,TMath::Sqrt(nofCINT7),
4984                                          nofCMUL7,TMath::Sqrt(nofCMUL7));
4985
4986     r->Set("Fnorm",nofCINT7/nofCINT7w0MUL,(nofCINT7/nofCINT7w0MUL)*AliAnalysisMuMuResult::ErrorAB( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
4987                                                                       nofCINT7,TMath::Sqrt(nofCINT7)));
4988
4989     Double_t yield =  r->GetValue("CorrNofJpsi") * mbeq;
4990
4991     Double_t yieldError = yield * AliAnalysisMuMuResult::ErrorAB( r->GetValue("CorrNofJpsi"), r->GetErrorStat("CorrNofJpsi"),
4992                                           mbeq,mbeqError);
4993
4994     r->Set(Form("Y%sJpsi",direction),yield,yieldError);
4995
4996     Double_t raa = yield/(tab*brsigmapp);
4997     Double_t raaError = AliAnalysisMuMuResult::ErrorABC(yield,yieldError,
4998                                                         tab,tabError,
4999                                                         brsigmapp,brsigmappError);
5000     r->Set(Form("R%sJpsi",direction),raa,raaError);
5001
5002     r->Set("NofInputJpsi",rsim->GetValue("NofInputJpsi",accEffSubResultName),rsim->GetErrorStat("NofInputJpsi",accEffSubResultName));
5003     r->Set("AccEffJpsi",rsim->GetValue("AccEffJpsi",accEffSubResultName),rsim->GetErrorStat("AccEffJpsi",accEffSubResultName));
5004
5005     AliAnalysisMuMuBinning::Range* bincm = new AliAnalysisMuMuBinning::Range(bin->What(),bin->Quantity(),ylowcms,yhighcms);
5006
5007     r->SetBin(*bincm);
5008
5009     finalBins.Add(bincm);
5010     finalResults.Add(r);
5011
5012     ++i;
5013   }
5014
5015   delete bins;
5016
5017   AliAnalysisMuMuSpectra* spectra = new AliAnalysisMuMuSpectra(type,direction);
5018
5019   for ( i = 0; i <= n; ++i )
5020   {
5021     Int_t j(i);
5022     if ( strcmp(direction,"pPb")==0 )
5023     {
5024       j = n-i;
5025     }
5026
5027     r = static_cast<AliAnalysisMuMuJpsiResult*>(finalResults.At(j));
5028
5029     bin = static_cast<AliAnalysisMuMuBinning::Range*>(finalBins.At(j));
5030
5031     spectra->AdoptResult(*bin,r);
5032   }
5033
5034
5035   delete corrSpectra;
5036
5037   return spectra;
5038 }
5039
5040 //_____________________________________________________________________________
5041 void AliAnalysisMuMu::SetConfig(const AliAnalysisMuMuConfig& config)
5042 {
5043   /// (re)set the config
5044   delete fConfig;
5045   fConfig = new AliAnalysisMuMuConfig(config);
5046 }
5047