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
  21 #include "AliAnalysisMuMuBinning.h"
  22 #include "AliAnalysisMuMuFnorm.h"
  23 #include "AliAnalysisMuMuGraphUtil.h"
  24 #include "AliAnalysisMuMuJpsiResult.h"
  25 #include "AliAnalysisMuMuSpectra.h"
  26 #include "AliAnalysisTriggerScalers.h"
  27 #include "AliCounterCollection.h"
  28 #include "AliHistogramCollection.h"
  29 #include "AliLog.h"
  30 #include "AliMergeableCollection.h"
  31 #include "Riostream.h"
  32 #include "TArrayL64.h"
  33 #include "TASImage.h"
  34 #include "TAxis.h"
  35 #include "TCanvas.h"
  36 #include "TColor.h"
  37 #include "TF1.h"
  38 #include "TFile.h"
  39 #include "TGraphErrors.h"
  40 #include "TGrid.h"
  41 #include "TH1.h"
  42 #include "TH2.h"
  43 #include "THashList.h"
  44 #include "TKey.h"
  45 #include "TLegend.h"
  46 #include "TLegendEntry.h"
  47 #include "TLine.h"
  48 #include "TList.h"
  49 #include "TMap.h"
  50 #include "TMath.h"
  51 #include "TObjArray.h"
  52 #include "TObjString.h"
  53 #include "TParameter.h"
  54 #include "TPaveText.h"
  55 #include "TROOT.h"
  56 #include "TStopwatch.h"
  57 #include "TStyle.h"
  58 #include "TSystem.h"
  59 #include <cassert>
  60 #include <map>
  61 #include <set>
  62 #include <string>
  63
  64 ClassImp(AliAnalysisMuMu)
  65
  66 TString AliAnalysisMuMu::fgOCDBPath("raw://");
  67
  68 TString AliAnalysisMuMu::fgDefaultDimuonTriggers("CMUL7-B-NOPF-MUON");
  69
  70 //,CMUL7-S-NOPF-ALLNOTRD,CMUL7-S-NOPF-MUON,CMUL8-S-NOPF-MUON,CMUL7-B-NOPF-ALLNOTRD,CMUU7-B-NOPF-ALLNOTRD,CMUU7-B-NOPF-MUON,CPBI1MUL-B-NOPF-MUON,CMULLO-B-NOPF-MUON");
  71
  72 TString AliAnalysisMuMu::fgDefaultMuonTriggers("CMSL7-S-NOPF-MUON");
  73
  74 //,CMSL7-S-NOPF-ALLNOTRD,CMSL8-S-NOPF-MUON,CMSL8-S-NOPF-ALLNOTRD,CMSL7-B-NOPF-MUON,CMUS1-B-NOPF-MUON,CMUS7-B-NOPF-MUON,CMSNGL-B-NOPF-MUON");
  75
  76 TString AliAnalysisMuMu::fgDefaultMinbiasTriggers("CINT7-B-NOPF-ALLNOTRD");
  77
  78 //,CINT7-S-NOPF-ALLNOTRD,CINT8-B-NOPF-ALLNOTRD,CINT8-S-NOPF-ALLNOTRD,CINT1-B-NOPF-ALLNOTRD,CPBI2_B1-B-NOPF-ALLNOTRD");
  79
  80 TString AliAnalysisMuMu::fgDefaultEventSelectionList("PSALL"); // for real data, for simulation see AliAnalysisMuMu ctor
  81
  82 TString AliAnalysisMuMu::fgDefaultPairSelectionList("pMATCHLOWRABSBOTH");
  83
  84 TString AliAnalysisMuMu::fgDefaultCentralitySelectionList("PP");
  85
  86 //TString AliAnalysisMuMu::fgDefaultFitTypeList("PSILOW:2,PSILOWalphaLow0.984nLow5.839alphaUp1.972nUp3.444:2,PSILOWMCTAILS:2");
  87 TString AliAnalysisMuMu::fgDefaultFitTypeList("PSILOWalphaLow0.984nLow5.839alphaUp1.972nUp3.444:2,PSILOWMCTAILS:2");
  88
  89 TString AliAnalysisMuMu::fgDefaultEventSelectionForSimulations("ALL");
  90 TString AliAnalysisMuMu::fgDefaultDimuonTriggerForSimulations("CMULLO-B-NOPF-MUON");
  91
  92 Bool_t AliAnalysisMuMu::fgIsCompactGraphs(kFALSE);
  93
  94 //_____________________________________________________________________________
  95 TString First(const TString s)
  96 {
  97   TString rv;
  98
  99   TObjArray* tokens = s.Tokenize(",");
 100
 101   if (!tokens) return rv;
 102
 103   rv = static_cast<TObjString*>(tokens->First())->String();
 104
 105   delete tokens;
 106
 107   return rv;
 108 }
 109
 110 //_____________________________________________________________________________
 111 TString FindTrigger(const AliMergeableCollection& mc,
 112                     const char* base,
 113                     const char* selection,
 114                     const char* paircut,
 115                     const char* centrality)
 116 {
 117   /// find the trigger containing the MinvPt histograms
 118
 119   std::vector<std::string> trigger2test;
 120
 121   //  trigger2test.push_back(Form("%s5-B-NOPF-ALLNOTRD",base));
 122   //  trigger2test.push_back(Form("%s1-B-NOPF-ALLNOTRD",base));
 123   //  trigger2test.push_back(Form("%s1B-ABCE-NOPF-MUON",base));
 124   if ( TString(base).Contains("||") || TString(base).Contains("-") )
 125   {
 126     trigger2test.push_back(base);
 127   }
 128   else
 129   {
 130     trigger2test.push_back(Form("%s-B-NOPF-ALLNOTRD",base));
 131     trigger2test.push_back(Form("%s-B-NOPF-MUON",base));
 132     trigger2test.push_back(Form("%s-S-NOPF-ALLNOTRD",base));
 133     trigger2test.push_back(Form("%s-S-NOPF-MUON",base));
 134   }
 135   trigger2test.push_back("ANY");
 136
 137   for ( std::vector<std::string>::size_type i = 0; i < trigger2test.size(); ++i )
 138   {
 139     std::string trigger = trigger2test[i];
 140
 141     if ( mc.GetObject(Form("/%s/%s/%s/%s",selection,trigger.c_str(),centrality,paircut),"MinvUS") ||
 142         mc.GetObject(Form("/%s/%s/%s/%s",selection,trigger.c_str(),centrality,paircut),"MinvUSPt")
 143         )
 144     {
 145       return trigger.c_str();
 146     }
 147   }
 148
 149   //  AliWarningGeneral("FindTrigger",Form("DID NOT FIND TRIGGER base=%s selection=%s paircut=%s centrality=%s",
 150   //                  base,selection,paircut,centrality));
 151   //  for ( std::vector<std::string>::size_type i = 0; i < trigger2test.size(); ++i )
 152   //  {
 153   //    AliWarningGeneral("FindTrigger",Form("tested trigger = %s",trigger2test[i].c_str()));
 154   //  }
 155   return "";
 156 }
 157
 158 //_____________________________________________________________________________
 159 AliAnalysisMuMu::AliAnalysisMuMu(const char* filename, const char* associatedSimFileName) : TObject(),
 160 fFilename(filename),
 161 fCounterCollection(0x0),
 162 fDimuonTriggers(fgDefaultDimuonTriggers),
 163 fMuonTriggers(fgDefaultMuonTriggers),
 164 fMinbiasTriggers(fgDefaultMinbiasTriggers),
 165 fEventSelectionList(fgDefaultEventSelectionList),
 166 fPairSelectionList(fgDefaultPairSelectionList),
 167 fCentralitySelectionList(fgDefaultCentralitySelectionList),
 168 fFitTypeList(fgDefaultFitTypeList),
 169 fBinning(0x0),
 170 fMergeableCollection(0x0),
 171 fRunNumbers(),
 172 fCorrectionPerRun(0x0),
 173 fAssociatedSimulation(0x0)
 174 {
 175   // ctor
 176
 177   GetCollections(fFilename,fMergeableCollection,fCounterCollection,fBinning,fRunNumbers);
 178
 179   if ( fCounterCollection )
 180   {
 181     if (IsSimulation())
 182     {
 183       SetEventSelectionList("ALL");
 184       SetDimuonTriggerList("CMULLO-B-NOPF-MUON");
 185       SetFitTypeList("COUNTJPSI:1");
 186 //    SetFitTypeList("PSI1:1,COUNTJPSI:1");
 187     }
 188
 189     if ( strlen(associatedSimFileName) )
 190     {
 191       fAssociatedSimulation = new AliAnalysisMuMu(associatedSimFileName);
 192     }
 193   }
 194 }
 195
 196 //_____________________________________________________________________________
 197 AliAnalysisMuMu::~AliAnalysisMuMu()
 198 {
 199   // dtor
 200
 201   if ( fAssociatedSimulation )
 202   {
 203     fAssociatedSimulation->Update();
 204   }
 205
 206   Update();
 207
 208   delete fCounterCollection;
 209   delete fBinning;
 210   delete fMergeableCollection;
 211   delete fCorrectionPerRun;
 212   delete fAssociatedSimulation;
 213 }
 214
 215 //_____________________________________________________________________________
 216 void AliAnalysisMuMu::BasicCounts(Bool_t detailTriggers,
 217                                   ULong64_t* totalNmb,
 218                                   ULong64_t* totalNmsl,
 219                                   ULong64_t* totalNmul)
 220 {
 221   // Report of some basic numbers, like number of MB and MUON triggers,
 222   // both before and after physics selection, and comparison with
 223   // the total number of such triggers (taken from the OCDB scalers)
 224   // if requested.
 225   //
 226   // filename is assumed to be a root filecontaining a list containing
 227   //    an AliCounterCollection (or directly an AliCounterCollection)
 228   //
 229   // if detailTriggers is kTRUE, each kind of (MB,MUL,MSL) is counted separately
 230   //
 231
 232   if (!fMergeableCollection || !fCounterCollection) return;
 233
 234   TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
 235   TIter nextRun(runs);
 236
 237   TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
 238   TIter nextTrigger(triggers);
 239
 240   TObjArray* events = fCounterCollection->GetKeyWords("event").Tokenize(",");
 241
 242   Bool_t doPS = (events->FindObject("PSALL") != 0x0);
 243
 244   TObjString* srun;
 245   TObjString* strigger;
 246
 247   ULong64_t localNmb(0);
 248   ULong64_t localNmsl(0);
 249   ULong64_t localNmul(0);
 250
 251   if ( totalNmb) *totalNmb = 0;
 252   if ( totalNmsl) *totalNmsl = 0;
 253   if ( totalNmul ) *totalNmul = 0;
 254
 255   while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
 256   {
 257     std::cout << Form("RUN %09d ",srun->String().Atoi());
 258
 259     TString details;
 260     ULong64_t nmb(0);
 261     ULong64_t nmsl(0);
 262     ULong64_t nmul(0);
 263
 264     nextTrigger.Reset();
 265
 266     Int_t nofPS(0);
 267
 268     while ( ( strigger = static_cast<TObjString*>(nextTrigger()) ) )
 269     {
 270
 271       if ( !fgDefaultMinbiasTriggers.Contains(strigger->String().Data()) &&
 272            !fgDefaultMuonTriggers.Contains(strigger->String().Data()) &&
 273            !fgDefaultDimuonTriggers.Contains(strigger->String().Data()) ) continue;
 274
 275       ULong64_t n = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
 276                                                     strigger->String().Data(),"ALL",srun->String().Atoi())));
 277
 278       details += TString::Format("\n%50s %10lld",strigger->String().Data(),n);
 279
 280
 281       ULong64_t nps = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
 282                                                       strigger->String().Data(),"PSALL",srun->String().Atoi())));
 283
 284       if ( doPS )
 285       {
 286         details += TString::Format(" PS %5.1f %%",nps*100.0/n);
 287       }
 288
 289       if (nps)
 290       {
 291         ++nofPS;
 292       }
 293
 294       if ( fMinbiasTriggers.Contains(strigger->String()) )
 295       {
 296         nmb += n;
 297         if ( totalNmb) (*totalNmb) += n;
 298         localNmb += n;
 299       }
 300       else if ( fMuonTriggers.Contains(strigger->String()) )
 301       {
 302         nmsl += n;
 303         if ( totalNmsl) (*totalNmsl) += n;
 304         localNmsl += n;
 305       }
 306       else if ( fDimuonTriggers.Contains(strigger->String()) )
 307       {
 308         nmul += n;
 309         if ( totalNmul ) (*totalNmul) += n;
 310         localNmul += n;
 311       }
 312     }
 313
 314     std::cout << Form("MB %10lld MSL %10lld MUL %10lld %s",
 315                  nmb,nmsl,nmul,(nofPS == 0 ? "(NO PS AVAIL)": ""));
 316
 317     if ( detailTriggers )
 318     {
 319       std::cout << details.Data();
 320     }
 321     std::cout << std::endl;
 322   }
 323
 324   if ( !totalNmul && !totalNmsl && !totalNmb )
 325   {
 326     std::cout << std::endl << Form("%13s MB %10lld MSL %10lld MUL %10lld ","TOTAL",
 327                                    localNmb,localNmsl,localNmul) << std::endl;
 328   }
 329
 330   delete runs;
 331   delete triggers;
 332   delete events;
 333 }
 334
 335 //_____________________________________________________________________________
 336 void AliAnalysisMuMu::BasicCountsEvolution(const char* filelist, Bool_t detailTriggers)
 337 {
 338   // Report of some basic numbers, like number of MB and MUON triggers,
 339   // both before and after physics selection, and comparison with
 340   // the total number of such triggers (taken from the OCDB scalers)
 341   // if requested.
 342   //
 343   // if detailTriggers is kTRUE, each kind of (MB,MUL,MSL) is counted separately
 344   //
 345   // To change the list of (single muon, dimuon, MB) triggers, use
 346   // the SetDefault*TriggerList methods prior to call this one
 347   //
 348
 349   TObjArray* files = ReadFileList(filelist);
 350
 351   if (!files || files->IsEmpty() ) return;
 352
 353   TIter next(files);
 354   TObjString* str;
 355
 356   ULong64_t totalNmb(0);
 357   ULong64_t totalNmsl(0);
 358   ULong64_t totalNmul(0);
 359
 360   while ( ( str = static_cast<TObjString*>(next()) ) )
 361   {
 362     AliAnalysisMuMu m(str->String().Data());
 363
 364     ULong64_t nmb(0);
 365     ULong64_t nmsl(0);
 366     ULong64_t nmul(0);
 367
 368     m.BasicCounts(detailTriggers,&nmb,&nmsl,&nmul);
 369
 370     totalNmb += nmb;
 371     totalNmsl += nmsl;
 372     totalNmul += nmul;
 373   }
 374
 375   std::cout << std::endl << Form("%13s MB %10lld MSL %10lld MUL %10lld ","TOTAL",
 376                                  totalNmb,totalNmsl,totalNmul) << std::endl;
 377
 378 }
 379
 380 //_____________________________________________________________________________
 381 void AliAnalysisMuMu::CentralityCheck(const char* filelist)
 382 {
 383   // Check if we get correctly filled centrality
 384
 385   TObjArray* files = ReadFileList(filelist);
 386
 387   if (!files || files->IsEmpty() ) return;
 388
 389   TIter next(files);
 390   TObjString* str;
 391
 392   while ( ( str = static_cast<TObjString*>(next()) ) )
 393   {
 394     AliMergeableCollection* mc(0x0);
 395     AliCounterCollection* cc(0x0);
 396     AliAnalysisMuMuBinning* bin(0x0);
 397     std::set<int> runnumbers;
 398
 399     if (!GetCollections(str->String().Data(),mc,cc,bin,runnumbers)) continue;
 400
 401     int run = RunNumberFromFileName(str->String().Data());
 402
 403     TH1* h = mc->Histo("/ALL/CPBI1MUL-B-NOPF-MUON/Centrality");
 404
 405     float percent(0);
 406
 407     if (h)
 408     {
 409       percent = 100*h->Integral(1,1)/h->Integral();
 410     }
 411
 412     std::cout << Form("RUN %09d PERCENT %7.2f",run,percent) << std::endl;
 413
 414     delete mc;
 415   }
 416
 417   gROOT->CloseFiles();
 418
 419 }
 420
 421 //_____________________________________________________________________________
 422 void AliAnalysisMuMu::CleanAllSpectra()
 423 {
 424   /// Delete all the spectra we may have
 425
 426   MC()->RemoveByType("AliAnalysisMuMuSpectra");
 427   Update();
 428 }
 429
 430
 431 //_____________________________________________________________________________
 432 TObjArray* AliAnalysisMuMu::CompareJpsiPerCMUUWithBackground(const char* jpsiresults,
 433                                                                    const char* backgroundresults)
 434 {
 435   TFile* fjpsi = FileOpen(jpsiresults);
 436   TFile* fbck = FileOpen(backgroundresults);
 437
 438   if (!fjpsi || !fbck) return 0x0;
 439
 440   TGraph* gjpsi = static_cast<TGraph*>(fjpsi->Get("jpsipercmuu"));
 441
 442   std::vector<std::string> checks;
 443
 444   checks.push_back("muminus-CMUU7-B-NOPF-ALLNOTRD");
 445   checks.push_back("muplus-CMUU7-B-NOPF-ALLNOTRD");
 446   checks.push_back("muminus-CMUSH7-B-NOPF-MUON");
 447   checks.push_back("muplus-CMUSH7-B-NOPF-MUON");
 448
 449   if (!gjpsi) return 0x0;
 450
 451   TObjArray* a = new TObjArray;
 452   a->SetOwner(kTRUE);
 453
 454   for ( std::vector<std::string>::size_type j = 0; j < checks.size(); ++j )
 455   {
 456
 457     TGraph* gback = static_cast<TGraph*>(fbck->Get(checks[j].c_str()));
 458
 459     if (!gback) continue;
 460
 461     if ( gjpsi->GetN() != gback->GetN() )
 462     {
 463       AliErrorClass("graphs have different number of points !");
 464       continue;
 465     }
 466
 467     TGraphErrors* g = new TGraphErrors(gjpsi->GetN());
 468
 469     for ( int i = 0; i < gjpsi->GetN(); ++i )
 470     {
 471       double r1,r2,y1,y2;
 472
 473       gjpsi->GetPoint(i,r1,y1);
 474       gback->GetPoint(i,r2,y2);
 475
 476       if ( r1 != r2 )
 477       {
 478         AliWarningClass(Form("run[%d]=%d vs %d",i,(int)r1,(int)r2));
 479         continue;
 480       }
 481
 482       g->SetPoint(i,y2,y1);
 483       //    g->SetPointError(i,gjpsi->GetErrorY(i),gback->GetErrorY(i));
 484     }
 485
 486     g->SetMarkerStyle(25+j);
 487     g->SetMarkerSize(1.2);
 488     if (j==0)
 489     {
 490       g->Draw("ap");
 491     }
 492     else
 493     {
 494       g->Draw("p");
 495     }
 496     g->SetLineColor(j+1);
 497     g->SetMarkerColor(j+1);
 498     g->SetName(checks[j].c_str());
 499     a->AddLast(g);
 500   }
 501
 502   return a;
 503 }
 504
 505 //_____________________________________________________________________________
 506 TGraph* AliAnalysisMuMu::CompareJpsiPerCMUUWithSimu(const char* realjpsiresults,
 507                                                              const char* simjpsiresults)
 508 {
 509   TFile* freal = FileOpen(realjpsiresults);
 510   TFile* fsim = FileOpen(simjpsiresults);
 511
 512   if (!freal || !fsim) return 0x0;
 513
 514   TGraph* greal = static_cast<TGraph*>(freal->Get("jpsipercmuu"));
 515   TGraph* gsim = static_cast<TGraph*>(fsim->Get("jpsipercmuu"));
 516
 517   TObjArray* a = new TObjArray;
 518   a->SetOwner(kTRUE);
 519
 520   if ( greal->GetN() != gsim->GetN() )
 521   {
 522     AliErrorClass("graphs have different number of points !");
 523     return 0x0;
 524   }
 525
 526   TGraphErrors* g = new TGraphErrors(greal->GetN());
 527   TGraphErrors* gratio = new TGraphErrors(greal->GetN());
 528
 529   for ( int i = 0; i < greal->GetN(); ++i )
 530   {
 531     double r1,r2,y1,y2;
 532
 533     greal->GetPoint(i,r1,y1);
 534     gsim->GetPoint(i,r2,y2);
 535
 536     if ( r1 != r2 )
 537     {
 538       AliWarningClass(Form("run[%d]=%d vs %d",i,(int)r1,(int)r2));
 539       continue;
 540     }
 541
 542     double ratio(0.0);
 543
 544     if ( TMath::Abs(y1)<1E-6 || TMath::Abs(y2)<1E-6)
 545     {
 546       g->SetPoint(i,0,0);
 547       g->SetPointError(i,0,0);
 548     }
 549     else
 550     {
 551       g->SetPoint(i,y2,y1);
 552       g->SetPointError(i,greal->GetErrorY(i),gsim ->GetErrorY(i));
 553       ratio = y2/y1;
 554     }
 555     gratio->SetPoint(i,r1,ratio);
 556   }
 557
 558   g->SetMarkerStyle(25);
 559   g->SetMarkerSize(1.2);
 560
 561   new TCanvas;
 562
 563   g->Draw("ap");
 564
 565   g->SetLineColor(1);
 566   g->SetMarkerColor(1);
 567   g->SetName("jpsipercmuurealvssim");
 568
 569   new TCanvas;
 570
 571   greal->Draw("alp");
 572   gsim->SetLineColor(4);
 573
 574   gsim->Draw("lp");
 575
 576   new TCanvas;
 577   gratio->Draw("alp");
 578
 579   return g;
 580 }
 581
 582 //_____________________________________________________________________________
 583 TObjArray* AliAnalysisMuMu::ComputeBackgroundEvolution(const char* filelist,
 584                                                        const char* triggerList,
 585                                                        Double_t ptmin,
 586                                                        const char* outputFile,
 587                                                        const char* outputMode)
 588 {
 589   // triggerList is a list of complete trigger names, separated by space
 590   // of the triggers to consider : only the first one found in the list
 591   // is used for each run.
 592
 593   TObjArray* files = ReadFileList(filelist);
 594
 595   if (!files || files->IsEmpty() ) return 0x0;
 596
 597   TIter next(files);
 598   TObjString* str;
 599
 600   const char* ps = "PSALL";
 601   const char* centrality = "PP";
 602   const char* ts1 = "sMATCHLOWRABS";
 603   const char* ts2 = "sMATCHLOWRABSDCA";
 604
 605   std::map<std::string, std::vector<float> > runs;
 606   std::map<std::string, std::vector<float> > errruns;
 607   std::map<std::string, std::vector<float> > yplus,erryplus;
 608   std::map<std::string, std::vector<float> > yminus,erryminus;
 609
 610   TObjArray* triggers = TString(triggerList).Tokenize(",");
 611
 612   TObjArray* a = new TObjArray;
 613   a->SetOwner(kTRUE);
 614
 615   Bool_t bothSigns(kFALSE);
 616
 617   while ( ( str = static_cast<TObjString*>(next()) ) )
 618   {
 619     AliInfoClass(str->String().Data());
 620
 621     AliMergeableCollection* mc(0x0);
 622     AliCounterCollection* cc(0x0);
 623     AliAnalysisMuMuBinning* bin(0x0);
 624     std::set<int> runnumbers;
 625
 626     if (!GetCollections(str->String().Data(),mc,cc,bin,runnumbers)) continue;
 627
 628     TIter nextObject(mc->CreateIterator());
 629     TObject* o;
 630     int nplus(0), nminus(0);
 631
 632     while ( (  o = nextObject() ) )
 633     {
 634       if ( o->InheritsFrom("TH1") )
 635       {
 636         continue;
 637       }
 638
 639       TH1* h = static_cast<TH1*>(o);
 640
 641       if ( TString(h->GetName()).EndsWith("Plus") )
 642       {
 643         nplus++;
 644       }
 645       if ( TString(h->GetName()).EndsWith("Minus") )
 646       {
 647         nminus++;
 648       }
 649     }
 650
 651     if  (nminus==nplus && nplus>0 )
 652     {
 653       bothSigns = kTRUE;
 654     }
 655
 656     AliInfoClass(Form("Both signs = %d",bothSigns));
 657
 658     TIter nextTrigger(triggers);
 659     TObjString* trigger;
 660     TH1* h1p(0x0);
 661     TH1* h1m(0x0);
 662     TH1* h2p(0x0);
 663     TH1* h2m(0x0);
 664     TString format;
 665
 666     while ( ( trigger = static_cast<TObjString*>(nextTrigger()) ) )
 667     {
 668       if  (bothSigns)
 669       {
 670         format = "/%s/%s/%s/%s/PtEtaMuPlus:py";
 671       }
 672       else
 673       {
 674         format = "/%s/%s/%s/%s/PtEtaMu:py";
 675       }
 676
 677       TString hname(Form(format.Data(),ps,trigger->String().Data(),centrality,ts1));
 678
 679       h1p = mc->Histo(hname.Data());
 680
 681       if (!h1p)
 682       {
 683         AliInfoClass(Form("Could not get %s",hname.Data()));
 684         continue;
 685       }
 686
 687       AliInfoClass(Form("Will use trigger %s",trigger->String().Data()));
 688
 689       h2p = mc->Histo(Form(format.Data(),ps,trigger->String().Data(),centrality,ts2));
 690
 691       if ( bothSigns )
 692       {
 693         format.ReplaceAll("Plus","Minus");
 694         h1m = mc->Histo(Form(format.Data(),ps,trigger->String().Data(),centrality,ts1));
 695         h2m = mc->Histo(Form(format.Data(),ps,trigger->String().Data(),centrality,ts2));
 696       }
 697       else
 698       {
 699         h2m=h2p;
 700         h1m=h1p;
 701       }
 702
 703       if (h1m && h2m && h1p && h2p)
 704       {
 705         Int_t bin1 = h1m->GetXaxis()->FindBin(ptmin);
 706         Int_t bin2 = h1m->GetXaxis()->GetNbins();
 707
 708         runs[trigger->String().Data()].push_back(RunNumberFromFileName(str->String().Data()));
 709         errruns[trigger->String().Data()].push_back(0.5);
 710
 711         double e1,e2;
 712         double v1 = h2m->IntegralAndError(bin1,bin2,e1);
 713         double v2 = h1m->IntegralAndError(bin1,bin2,e2);
 714         double value = 100*(1.0-v1/v2);
 715         e1/=v1;
 716         e2/=v2;
 717         yminus[trigger->String().Data()].push_back(value);
 718 //        double e1 = 1.0/TMath::Sqrt(h1m->GetEntries());
 719 //        double e2 = 1.0/TMath::Sqrt(h2m->GetEntries());
 720         erryminus[trigger->String().Data()].push_back(TMath::Sqrt(e1*e1+e2*e2)*value);
 721
 722         v1=h2p->IntegralAndError(bin1,bin2,e1);
 723         v2=h1p->IntegralAndError(bin1,bin2,e1);
 724         value = 100*(1.0-v1/v2);
 725         e1/=v1;
 726         e2/=v2;
 727         yplus[trigger->String().Data()].push_back(value);
 728 //        e1 = 1.0/TMath::Sqrt(h1p->GetEntries());
 729 //        e2 = 1.0/TMath::Sqrt(h2p->GetEntries());
 730         erryplus[trigger->String().Data()].push_back(TMath::Sqrt(e1*e1+e2*e2)*value);
 731       }
 732       else
 733       {
 734         std::cout << Form("Error : h1m %p h2m %p h1p %p h2p %p",h1m,h2m,h1p,h2p) << std::endl;
 735       }
 736     }
 737
 738     delete mc;
 739     delete cc;
 740     TFile* f = static_cast<TFile*>(gROOT->GetListOfFiles()->FindObject(str->String().Data()));
 741     delete f;
 742   }
 743
 744   delete triggers;
 745
 746   TFile* f = new TFile(outputFile,outputMode);
 747
 748   std::map<std::string, std::vector<float> >::const_iterator it;
 749
 750   for ( it = runs.begin(); it != runs.end(); ++it )
 751   {
 752     std::string triggerName = it->first;
 753
 754   TGraphErrors* gp = new TGraphErrors(runs[triggerName].size(),&runs[triggerName][0],&yplus[triggerName][0],&errruns[triggerName][0],&erryplus[triggerName][0]);
 755   TGraphErrors* gm(0x0);
 756
 757   if ( bothSigns )
 758   {
 759     gm = new TGraphErrors(runs[triggerName].size(),&runs[triggerName][0],&yminus[triggerName][0],&errruns[triggerName][0],&erryminus[triggerName][0]);
 760   }
 761
 762   if ( bothSigns )
 763   {
 764     gp->Write(Form("muplus_%s",triggerName.c_str()),TObject::kOverwrite);
 765     gm->Write(Form("muminus_%s",triggerName.c_str()),TObject::kOverwrite);
 766   }
 767   else
 768   {
 769     gp->Write(Form("mu_%s",triggerName.c_str()),TObject::kOverwrite);
 770   }
 771
 772   }
 773
 774   delete f;
 775
 776   return a;
 777 }
 778
 779 //_____________________________________________________________________________
 780 TMap*
 781 AliAnalysisMuMu::ComputeJpsiEvolution(const char* filelist, const char* triggerList,
 782                                       const char* outputFile)
 783 {
 784   /// Compute some jpsi information for a list of files / trigger combinations
 785
 786   TObjArray* files = ReadFileList(filelist);
 787
 788   if (!files || files->IsEmpty() ) return 0x0;
 789
 790   TMap results; // one TObjString->TObjArray per file
 791   results.SetOwnerKeyValue(kTRUE,kTRUE);
 792
 793   TIter nextFile(files);
 794   TObjString* str;
 795   TString fitType;
 796
 797 //  while ( ( str = static_cast<TObjString*>(nextFile()) ) )
 798 //  {
 799 //    std::cout << str->String().Data() << std::endl;
 800 //
 801 //    AliAnalysisMuMu m(str->String().Data());
 802 //
 803 //    m.SetDimuonTriggerList(triggerList);
 804 //
 805 //    TMap* map = m.Jpsi();
 806 //
 807 //    if (!map)
 808 //    {
 809 //      AliWarningClass(Form("Got no jpsi for %s",str->String().Data()));
 810 //    }
 811 //
 812 //    results.Add(new TObjString(str->String()), map);
 813 //  }
 814 //
 815 //  if (!results.GetSize()) return 0x0;
 816
 817   // compute the total over all files
 818
 819   TMap* total = new TMap;
 820   total->SetOwnerKeyValue(kTRUE,kTRUE);
 821
 822   nextFile.Reset();
 823   TObjArray* triggers = TString(triggerList).Tokenize(",");
 824
 825   TIter nextTrigger(triggers);
 826   TObjString* trigger(0x0);
 827
 828   while ( ( trigger = static_cast<TObjString*>(nextTrigger())))
 829   {
 830     nextFile.Reset();
 831
 832     TList l;
 833     AliAnalysisMuMuResult* ref(0x0);
 834
 835     while ( ( str = static_cast<TObjString*>(nextFile()) ) )
 836     {
 837 //      TObjArray* a = static_cast<TObjArray*>(results.GetValue(str->String().Data()));
 838
 839       AliInfoClass("FIXME: write the merging of AliAnalysisMuMuResult objects !");
 840
 841 //      AliAnalysisMuMuResult* r(0x0);
 842 //
 843 //      if (a)
 844 //      {
 845 //        r = static_cast<AliAnalysisMuMuResult*>(a->FindObject(trigger->String().Data()));
 846 //
 847 //        if (r)
 848 //        {
 849 //          if (!ref) ref = r;
 850 //
 851 //          if ( !hminv )
 852 //          {
 853 //            TH1* htmp = static_cast<TH1*>(r->Minv()->Clone());
 854 //            if (!htmp)
 855 //            {
 856 //              continue;
 857 //            }
 858 //            hminv = htmp;
 859 //          }
 860 //          else
 861 //          {
 862 //            l.Add(r->Minv());
 863 //          }
 864 //
 865 //          n += r->NofTriggers();
 866 //          ++nruns;
 867 //        }
 868 //      }
 869     }
 870
 871 //    sum->Merge(&l);
 872
 873     if (!ref) continue;
 874
 875
 876 //    AliAnalysisMuMuResult* sum = new AliAnalysisMuMuResult(*hminv,
 877 //                                                           ref->TriggerClass(),
 878 //                                                           ref->EventSelection(),
 879 //                                                           ref->PairSelection(),
 880 //                                                           ref->CentralitySelection(),
 881 //                                                           AliAnalysisMuMuBinning::Range());
 882 //
 883 //    sum->SetNofTriggers(n);
 884 //
 885 //    sum->SetNofRuns(nruns);
 886 //
 887 //    sum->Fit(1);
 888 //
 889 //    total->Add(new TObjString(trigger->String().Data()),sum);
 890
 891   }
 892
 893   AliInfoClass("--------------------------------------");
 894   StdoutToAliInfoClass(total->Print(););
 895
 896   AliInfoClass("---------------------------Going to write file");
 897
 898   TFile* fout = new TFile(outputFile,"RECREATE");
 899
 900   results.Write("rbr",TObject::kSingleKey);
 901
 902   total->Write("total",TObject::kSingleKey);
 903
 904   fout->Close();
 905
 906   delete fout;
 907
 908   AliInfoClass(Form("%d files analyzed",files->GetEntries()));
 909
 910   return total;
 911 }
 912
 913 //_____________________________________________________________________________
 914 Bool_t AliAnalysisMuMu::DecodeFileName(const char* filename,
 915                                              TString& period,
 916                                              int& esdpass,
 917                                              int& aodtrain,
 918                                              int& runnumber)
 919 {
 920   // tries to extract period and pass numbers from a file name.
 921
 922   esdpass=aodtrain=runnumber=-1;
 923   period="";
 924
 925   TString sfile(gSystem->BaseName(filename));
 926
 927   if (!sfile.BeginsWith("LHC") && !sfile.BeginsWith("SIM") )
 928   {
 929     // does not look nice but let's try at least to get the period
 930
 931     TObjArray* tmp = sfile.Tokenize(".");
 932
 933     for ( Int_t j = 0; j < tmp->GetEntries(); ++j )
 934     {
 935       TString s = static_cast<TObjString*>(tmp->At(j))->String();
 936       s.ToLower();
 937       if ( s.BeginsWith("lhc") )
 938       {
 939         period = s;
 940         period.ReplaceAll("lhc","LHC");
 941         break;
 942       }
 943     }
 944
 945     delete tmp;
 946
 947     if ( period.Length() )
 948     {
 949       return kTRUE;
 950     }
 951
 952     std::cerr << Form("filename %s does not start with LHC or SIM",filename) << std::endl;
 953     return kFALSE;
 954   }
 955
 956   int year;
 957   char p;
 958   Bool_t ok(kFALSE);
 959
 960   if ( sfile.BeginsWith("LHC") )
 961   {
 962     if ( sfile.Contains("pass") && sfile.Contains("AODMUON") )
 963     {
 964       int pass;
 965       sscanf(sfile.Data(),"LHC%2d%c_pass%d_AODMUON%03d_%09d",&year,&p,&pass,&aodtrain,&runnumber);
 966       period = TString::Format("LHC%2d%c",year,p);
 967       ok=kTRUE;
 968     }
 969     else if ( sfile.Contains("pass") && sfile.Contains("_muon_") && sfile.Contains("AOD000") )
 970     {
 971       // LHC11c_pass2_muon_AOD000_000152087.saf.root
 972       sscanf(sfile.Data(),"LHC%2d%c_pass%d_muon_AOD000_%09d",&year,&p,&esdpass,&runnumber);
 973       period = TString::Format("LHC%2d%c",year,p);
 974       ok=kTRUE;
 975     }
 976     else if ( sfile.Contains("_muon_calo_") && sfile.Contains("AODMUON000") )
 977     {
 978       //      LHC12h_muon_calo_AODMUON000_000190112.saf.root
 979       sscanf(sfile.Data(),"LHC%2d%c_muon_calo_AODMUON000_%09d",&year,&p,&runnumber);
 980       period = TString::Format("LHC%2d%c",year,p);
 981       ok=kTRUE;
 982     }
 983     else if ( sfile.Contains("_muon_calo_") && sfile.Contains("AOD000") )
 984     {
 985       //      LHC12h_muon_calo_AOD000_000190112.saf.root
 986       sscanf(sfile.Data(),"LHC%2d%c_muon_calo_AOD000_%09d",&year,&p,&runnumber);
 987       period = TString::Format("LHC%2d%c",year,p);
 988       ok=kTRUE;
 989     }
 990     else if ( sfile.Contains("AODMUON" ) )
 991     {
 992       sscanf(sfile.Data(),"LHC%2d%c_AODMUON%03d_%09d",&year,&p,&aodtrain,&runnumber);
 993       period = TString::Format("LHC%2d%c",year,p);
 994       ok=kTRUE;
 995     }
 996     else if ( sfile.Contains("AOD000") )
 997     {
 998       sscanf(sfile.Data(),"LHC%2d%c_muons_AOD000_%09d",&year,&p,&runnumber);
 999       period = TString::Format("LHC%2d%c",year,p);
1000       ok=kTRUE;
1001     }
1002     else if ( sfile.Contains("ESD_OUTER000"))
1003     {
1004       sscanf(sfile.Data(),"LHC%2d%c_cpass1_ESD_OUTER000_%09d",&year,&p,&runnumber);
1005       ok=kTRUE;
1006     }
1007   }
1008   else if ( sfile.BeginsWith("SIM_JPSI3" ) )
1009   {
1010     sscanf(sfile.Data(),"SIM_JPSI3_%09d",&runnumber);
1011     ok = kTRUE;
1012   }
1013   else if ( sfile.BeginsWith("SIM_UPSILON" ) )
1014   {
1015     sscanf(sfile.Data(),"SIM_UPSILON_%09d",&runnumber);
1016     ok = kTRUE;
1017   }
1018   else if ( sfile.BeginsWith("SIM_JPSI" ) )
1019   {
1020     sscanf(sfile.Data(),"SIM_JPSI_LHC%2d%c_%09d",&year,&p,&runnumber);
1021     period = TString::Format("LHC%2d%c",year,p);
1022     ok = kTRUE;
1023   }
1024
1025   if (!ok)
1026   {
1027     std::cerr << Form("Can not decode %s",filename) << std::endl;
1028     return kFALSE;
1029   }
1030
1031   return kTRUE;
1032 }
1033
1034 //_____________________________________________________________________________
1035 void AliAnalysisMuMu::DrawMinv(const char* type,
1036                                const char* particle,
1037                                const char* trigger,
1038                                const char* eventType,
1039                                const char* pairCut,
1040                                const char* centrality,
1041                                const char* subresultname,
1042                                const char* flavour) const
1043 {
1044   /// Draw minv spectra for binning of given type
1045
1046   if (!MC() || !BIN()) return;
1047
1048   TObjArray* bins = BIN()->CreateBinObjArray(particle,type,flavour);
1049   if (!bins)
1050   {
1051     AliError(Form("Could not get %s bins",type));
1052     return;
1053   }
1054
1055   Double_t xmin(-1);
1056   Double_t xmax(-1);
1057
1058   TString sparticle(particle);
1059   if ( sparticle=="PSI" )
1060   {
1061     xmin = 2;
1062     xmax = 6;
1063   }
1064
1065   Int_t nx(1);
1066   Int_t ny(1);
1067
1068   Int_t n = bins->GetEntries();
1069
1070   if ( n == 2 )
1071   {
1072     nx = 2;
1073   }
1074   else if ( n > 2 )
1075   {
1076     ny = TMath::Nint(TMath::Sqrt(n));
1077     nx = n/ny;
1078   }
1079
1080   TString stype(type);
1081   stype.ToUpper();
1082
1083   TString spectraName(Form("/%s/%s/%s/%s/%s-%s",eventType,trigger,centrality,pairCut,particle,stype.Data()));
1084
1085   if ( strlen(flavour))
1086   {
1087     spectraName += "-";
1088     spectraName += flavour;
1089   }
1090
1091   AliAnalysisMuMuSpectra* spectra = static_cast<AliAnalysisMuMuSpectra*>(MC()->GetObject(spectraName.Data()));
1092
1093   AliDebug(1,Form("spectraName=%s spectra=%p",spectraName.Data(),spectra));
1094
1095   TObjArray* spectraBins(0x0);
1096   if ( spectra )
1097   {
1098     spectraBins = spectra->BinContentArray();
1099   }
1100
1101   TCanvas* c = new TCanvas;
1102   c->Divide(nx,ny);
1103   c->Draw();
1104   gStyle->SetOptFit(1112);
1105
1106   c->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", "AliAnalysisMuMu",
1107               (void*)this, "ExecuteCanvasEvent(Int_t,Int_t,Int_t,TObject*)");
1108
1109
1110   TIter next(bins);
1111   AliAnalysisMuMuBinning::Range* r;
1112   Int_t ci(0);
1113
1114   while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(next())) )
1115   {
1116     TString name(Form("/%s/%s/%s/%s/MinvUS%s",eventType,trigger,centrality,pairCut,r->AsString().Data()));
1117
1118     AliDebug(1,name.Data());
1119
1120     AliAnalysisMuMuJpsiResult* spectraBin(0x0);
1121
1122     if ( spectraBins )
1123     {
1124       AliAnalysisMuMuResult* sr = static_cast<AliAnalysisMuMuResult*>(spectraBins->At(ci));
1125
1126       spectraBin = static_cast<AliAnalysisMuMuJpsiResult*>(sr->SubResult(subresultname));
1127
1128       AliDebug(1,Form("spectraBin(%s)=%p",subresultname,spectraBin));
1129     }
1130
1131     TH1* h = MC()->Histo(name.Data());
1132
1133     if ( spectraBin )
1134     {
1135       h = spectraBin->Minv();
1136     }
1137
1138     if (h)
1139     {
1140       ++ci;
1141       c->cd(ci);
1142       gPad->SetLogy();
1143       if (xmin>0)
1144       {
1145         h->GetXaxis()->SetRangeUser(xmin,xmax);
1146       }
1147       h->Draw("histes");
1148
1149       TObject* f1 = h->GetListOfFunctions()->FindObject("fitTotal");
1150       if (f1)
1151       {
1152         f1->Draw("same");
1153       }
1154
1155       gPad->Modified();
1156       gPad->Update();
1157
1158       TObject* stats = h->FindObject("stats");
1159       if (stats)
1160       {
1161         stats->Draw("same");
1162       }
1163     }
1164   }
1165
1166   delete bins;
1167 }
1168
1169 //_____________________________________________________________________________
1170 void AliAnalysisMuMu::DrawMinv(const char* type, const char* particle, const char* flavour, const char* subresultname) const
1171 {
1172   /// Draw minv spectra for binning of given type
1173
1174   DrawMinv(type,particle,
1175            First(DimuonTriggerList()).Data(),
1176            First(EventSelectionList()).Data(),
1177            First(PairSelectionList()).Data(),
1178            First(CentralitySelectionList()).Data(),
1179            subresultname,
1180            flavour);
1181 }
1182
1183 //___________________________________________________________________
1184 void AliAnalysisMuMu::ExecuteCanvasEvent(Int_t event, Int_t /*px*/, Int_t /*py*/, TObject *sel)
1185 {
1186   // Actions in reponse to mouse button events.
1187
1188   TCanvas* c = static_cast<TCanvas*>(gTQSender);
1189   TPad* pad = static_cast<TPad*>(c->GetSelectedPad());
1190   if (!pad) return;
1191
1192 //  if ((event == kButton1Down) ||
1193   if (event == kButton1Double)
1194   {
1195
1196 //    Float_t x = pad->AbsPixeltoX(px);
1197 //    Float_t y = pad->AbsPixeltoY(py);
1198 //    x = pad->PadtoX(x);
1199 //    y = pad->PadtoY(y);
1200
1201 //    std::cout << "event=" << event << " px=" << px << " py=" << py << " ";
1202
1203     if ( sel && sel->InheritsFrom("TH1") )
1204     {
1205       TCanvas* clocal = new TCanvas;
1206       clocal->SetLogy();
1207       clocal->Draw();
1208       sel->Draw();
1209     }
1210     else
1211     {
1212       TList* list = pad->GetListOfPrimitives();
1213       TIter next(list);
1214       TObject* h;
1215
1216       while ( ( h = next() ) )
1217       {
1218         if ( h->InheritsFrom("TH1") )
1219         {
1220           TCanvas* clocal = new TCanvas;
1221           clocal->SetLogy();
1222           clocal->Draw();
1223           h->Draw();
1224           break;
1225         }
1226       }
1227
1228     }
1229
1230 //      std::cout  << std::endl;
1231
1232       pad->Modified();
1233   }
1234
1235 }
1236
1237 //_____________________________________________________________________________
1238 TString
1239 AliAnalysisMuMu::ExpandPathName(const char* file)
1240 {
1241   // An expand method that lives alien URL as they are
1242   TString sfile;
1243
1244   if ( !sfile.BeginsWith("alien://") )
1245   {
1246     return gSystem->ExpandPathName(file);
1247   }
1248   else
1249   {
1250     if (!gGrid) TGrid::Connect("alien://");
1251     if (!gGrid) return "";
1252   }
1253
1254   return file;
1255 }
1256
1257 //_____________________________________________________________________________
1258 void AliAnalysisMuMu::TwikiOutputFnorm(const char* series) const
1259 {
1260   // make a twiki-compatible output of the Fnorm factor(s)
1261   TObjArray* what = TString(series).Tokenize(",");
1262   TObjString* s;
1263   TObjArray graphs;
1264   TIter next(what);
1265
1266   std::cout << "| *Run* |";
1267   while ( ( s = static_cast<TObjString*>(next())) )
1268   {
1269     TGraph* g = static_cast<TGraph*>(MC()->GetObject(Form("/FNORM/GRAPHS/%s",s->String().Data())));
1270     if (!g)
1271     {
1272       AliError(Form("Could not find graph for %s",s->String().Data()));
1273       continue;
1274     }
1275     std::cout << " *" << s->String().Data();
1276     if ( s->String().BeginsWith("RelDif") ) std::cout << " %";
1277     std::cout << "*|";
1278     graphs.Add(g);
1279   }
1280
1281   std::cout << std::endl;
1282
1283   TGraphErrors* g0 = static_cast<TGraphErrors*>(graphs.First());
1284   if (!g0) return;
1285
1286   for ( Int_t i = 0; i < g0->GetN(); ++i )
1287   {
1288     TString msg;
1289
1290     msg.Form("|%6d|",TMath::Nint(g0->GetX()[i]));
1291
1292     for ( Int_t j = 0; j < graphs.GetEntries(); ++j )
1293     {
1294       TGraphErrors* g = static_cast<TGraphErrors*>(graphs.At(j));
1295
1296       msg += TString::Format(" %6.2f +- %6.2f |",g->GetY()[i],g->GetEY()[i]);
1297     }
1298
1299     std::cout << msg.Data() << std::endl;
1300   }
1301
1302   next.Reset();
1303
1304   std::cout << "|*Weigthed mean (*)*|";
1305
1306   AliAnalysisMuMuResult* r = static_cast<AliAnalysisMuMuResult*>(MC()->GetObject("/FNORM/RESULTS/Fnorm"));
1307
1308   if (!r)
1309   {
1310     AliError("Could not find Fnorm result !");
1311     return;
1312   }
1313
1314
1315   while ( ( s = static_cast<TObjString*>(next())) )
1316   {
1317     TString var("Fnorm");
1318     TString unit;
1319
1320     if ( s->String().BeginsWith("Fnorm") )
1321     {
1322       r = static_cast<AliAnalysisMuMuResult*>(MC()->GetObject("/FNORM/RESULTS/Fnorm"));
1323     }
1324     else if ( s->String().BeginsWith("RelDif") )
1325     {
1326       r = static_cast<AliAnalysisMuMuResult*>(MC()->GetObject("/FNORM/RESULTS/RelDif"));
1327       unit = "%";
1328     }
1329
1330     r->Exclude("*");
1331     r->Include(s->String().Data());
1332
1333     std::cout << Form(" * %5.2f +- %5.2f %s * |",
1334                       r->GetValue(var.Data()),
1335                       r->GetErrorStat(var.Data()),
1336                       unit.Data());
1337   }
1338
1339   next.Reset();
1340
1341   std::cout << std::endl;
1342
1343   std::cout << "|*RMS*|";
1344
1345   while ( ( s = static_cast<TObjString*>(next())) )
1346   {
1347     TString var("Fnorm");
1348
1349     if ( s->String().BeginsWith("Fnorm") )
1350     {
1351       r = static_cast<AliAnalysisMuMuResult*>(MC()->GetObject("/FNORM/RESULTS/Fnorm"));
1352     }
1353     else if ( s->String().BeginsWith("RelDif") )
1354     {
1355       r = static_cast<AliAnalysisMuMuResult*>(MC()->GetObject("/FNORM/RESULTS/RelDif"));
1356     }
1357
1358     r->Exclude("*");
1359     r->Include(s->String().Data());
1360
1361     Double_t d = 100.0*r->GetRMS(var.Data())/r->GetValue(var.Data());
1362
1363     std::cout << Form(" * %5.2f (%5.2f %%) * |",
1364                       r->GetRMS(var.Data()),d);
1365   }
1366
1367   std::cout << std::endl;
1368   std::cout << "(*) weight is the number of CMUL7-B-NOPF-MUON triggers (physics-selected and pile-up corrected) in each run" << std::endl;
1369
1370   delete what;
1371 }
1372
1373 //_____________________________________________________________________________
1374 void AliAnalysisMuMu::FigureOutputFnorm(const char* filelist)
1375 {
1376   /// Make some figure of the Fnorm factors for files in filelist
1377
1378   TObjArray* periods = ReadFileList(filelist);
1379
1380   if (!periods || periods->IsEmpty() ) return;
1381
1382   TIter next(periods);
1383
1384   TObjArray fnormoffline1;
1385   TObjArray fnormoffline2;
1386   TObjArray reldif;
1387   TObjArray correctionPSMUL;
1388   TObjArray correctionPSMB;
1389   TObjArray correctionPUPS;
1390   TObjArray correctionPSRatio;
1391
1392   fnormoffline1.SetOwner(kTRUE);
1393   fnormoffline2.SetOwner(kTRUE);
1394   reldif.SetOwner(kTRUE);
1395   correctionPUPS.SetOwner(kTRUE);
1396   correctionPSMUL.SetOwner(kTRUE);
1397   correctionPSMB.SetOwner(kTRUE);
1398   correctionPSRatio.SetOwner(kTRUE);
1399
1400   for ( Int_t i = 0; i <= periods->GetLast(); ++i )
1401   {
1402     TString period("unknown");
1403
1404     TString filename(static_cast<TObjString*>(periods->At(i))->String());
1405
1406     Int_t dummy(0);
1407
1408     if (!DecodeFileName(filename,period,dummy,dummy,dummy))
1409     {
1410       continue;
1411     }
1412
1413     if ( gSystem->AccessPathName(filename) )
1414     {
1415       AliErrorClass(Form("Could not find file %s. Skipping it.",filename.Data()));
1416       continue;
1417
1418     }
1419
1420     AliAnalysisMuMu m(filename.Data());
1421
1422     fnormoffline1.Add(m.MC()->GetObject("/FNORM/GRAPHS/FnormOffline1PUPS")->Clone());
1423     fnormoffline2.Add(m.MC()->GetObject("/FNORM/GRAPHS/FnormOffline2PUPS")->Clone());
1424
1425
1426     correctionPSMUL.Add(m.MC()->GetObject("/FNORM/GRAPHS/CorrectionPSMUL")->Clone());
1427     correctionPSMB.Add(m.MC()->GetObject("/FNORM/GRAPHS/CorrectionPSMB")->Clone());
1428     correctionPUPS.Add(m.MC()->GetObject("/FNORM/GRAPHS/CorrectionPUPSMB")->Clone());
1429
1430     correctionPSRatio.Add(m.MC()->GetObject("/FNORM/GRAPHS/CorrectionPSRatio")->Clone());
1431
1432
1433     reldif.Add(m.MC()->GetObject("/FNORM/GRAPHS/RelDifFnormScalersPUPSvsFnormOffline2PUPS")->Clone());
1434
1435   }
1436
1437   AliAnalysisMuMuGraphUtil gu(fgOCDBPath);
1438
1439   gu.ShouldDrawPeriods(kTRUE);
1440
1441   TObjArray a;
1442
1443   a.Add(gu.Combine(fnormoffline1,fgIsCompactGraphs));
1444   a.Add(gu.Combine(fnormoffline2,fgIsCompactGraphs));
1445
1446   Double_t ymin(0.0);
1447   Double_t ymax(3000.0);
1448
1449   new TCanvas("fnormoffline","fnormoffline");
1450
1451   gu.PlotSameWithLegend(a,ymin,ymax);
1452
1453   new TCanvas("corrections","corrections");
1454
1455   a.Clear();
1456
1457   a.Add(gu.Combine(correctionPSMB,fgIsCompactGraphs));
1458   a.Add(gu.Combine(correctionPSMUL,fgIsCompactGraphs));
1459   a.Add(gu.Combine(correctionPUPS,fgIsCompactGraphs));
1460
1461   gu.PlotSameWithLegend(a,0.8,1.2);
1462
1463   new TCanvas("psratio","psratio");
1464
1465   a.Clear();
1466
1467   a.Add(gu.Combine(correctionPSRatio,fgIsCompactGraphs));
1468
1469   gu.PlotSameWithLegend(a,0.8,1.2);
1470
1471   new TCanvas("offlinevsscalers","fig:offlinevsscalers");
1472
1473   a.Clear();
1474
1475   a.Add(gu.Combine(reldif,fgIsCompactGraphs));
1476
1477   gu.PlotSameWithLegend(a,-15,15);
1478 }
1479
1480 //_____________________________________________________________________________
1481 TFile*
1482 AliAnalysisMuMu::FileOpen(const char* file)
1483 {
1484   // Open a file after expansion of its name
1485
1486   return TFile::Open(ExpandPathName(file).Data());
1487 }
1488
1489 //_____________________________________________________________________________
1490 TString AliAnalysisMuMu::First(const TString& list) const
1491 {
1492   TObjArray* a = list.Tokenize(",");
1493   if ( a->GetLast() < 0 ) return "";
1494
1495   TString rv = static_cast<TObjString*>(a->First())->String();
1496
1497   delete a;
1498
1499   return rv;
1500 }
1501
1502 //_____________________________________________________________________________
1503 AliAnalysisMuMuSpectra*
1504 AliAnalysisMuMu::FitParticle(const char* particle,
1505                              const char* trigger,
1506                              const char* eventType,
1507                              const char* pairCut,
1508                              const char* centrality,
1509                              const AliAnalysisMuMuBinning& binning)
1510 {
1511   // Fit the minv spectra to find the given particle
1512   // Returns an array of AliAnalysisMuMuResult objects
1513
1514   static int n(0);
1515
1516   TObjArray* bins = binning.CreateBinObjArray(particle);
1517   if (!bins)
1518   {
1519     AliError(Form("Did not get any bin for particle %s",particle));
1520     return 0x0;
1521   }
1522
1523   TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
1524   if ( !triggers->FindObject(trigger) )
1525   {
1526     AliDebug(1,Form("Did not find trigger %s",trigger));
1527     delete bins;
1528     delete triggers;
1529     return 0x0;
1530   }
1531   delete triggers;
1532
1533   TObjArray* events = fCounterCollection->GetKeyWords("event").Tokenize(",");
1534   if ( !events->FindObject(eventType) )
1535   {
1536     AliError(Form("Did not find eventType %s",eventType));
1537     delete bins;
1538     delete events;
1539     return 0x0;
1540   }
1541   delete events;
1542
1543   Int_t ntrigger = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s",trigger,eventType)));
1544
1545   if  (ntrigger<=0)
1546   {
1547     AliError(Form("No trigger for trigger:%s/event:%s",trigger,eventType));
1548     delete bins;
1549     return 0x0;
1550   }
1551
1552   TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
1553   Int_t nruns = runs->GetEntries();
1554   delete runs;
1555
1556
1557 //  binning.Print();
1558
1559   AliAnalysisMuMuSpectra* spectra(0x0);
1560
1561   AliAnalysisMuMuBinning::Range* bin;
1562   TIter next(bins);
1563
1564   TObjArray* fitTypeArray = fFitTypeList.Tokenize(",");
1565   TIter nextFitType(fitTypeArray);
1566   TObjString* fitType;
1567   TString flavour;
1568
1569   while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(next())) )
1570   {
1571     TString hname(Form("MinvUS%s",bin->AsString().Data()));
1572
1573     TH1* hminv = fMergeableCollection->Histo(Form("/%s/%s/%s/%s",eventType,trigger,centrality,pairCut),hname.Data());
1574
1575     if (!hminv)
1576     {
1577       if (!fBinning && bin->IsIntegrated() )
1578       {
1579         // old file, we only had MinvUSPt
1580         hminv = fMergeableCollection->Histo(Form("/%s/%s/%s/%s",eventType,trigger,centrality,pairCut),"MinvUSPt:py");
1581       }
1582
1583       if (!hminv)
1584       {
1585         AliDebug(1,Form("Could not find histo %s",hname.Data()));
1586         continue;
1587       }
1588     }
1589
1590     hminv = static_cast<TH1*>(hminv->Clone(Form("minv%d",n++)));
1591
1592
1593     AliAnalysisMuMuJpsiResult* r = new AliAnalysisMuMuJpsiResult(*hminv,
1594                                                                  trigger,
1595                                                                  eventType,
1596                                                                  pairCut,
1597                                                                  centrality,
1598                                                                  *bin);
1599
1600     r->SetNofTriggers(ntrigger);
1601     r->SetNofRuns(nruns);
1602
1603     nextFitType.Reset();
1604
1605     while ( ( fitType = static_cast<TObjString*>(nextFitType())) )
1606     {
1607       AliDebug(1,Form("<<<<<< fitType=%s bin=%s",fitType->String().Data(),bin->Flavour().Data()));
1608
1609       if ( fitType->String().BeginsWith("PSILOWMCTAILS") )
1610       {
1611         std::vector<Double_t> par;
1612         par = GetMCCB2Tails(*bin);
1613         if (!par.empty())
1614         {
1615           r->AddFit(fitType->String().Data(),par.size(),&par[0]);
1616         }
1617       }
1618       else
1619       {
1620         r->AddFit(fitType->String());
1621       }
1622     }
1623
1624     flavour = bin->Flavour();
1625
1626     if (!spectra)
1627     {
1628       TString spectraName(binning.GetName());
1629       if ( flavour.Length() > 0 )
1630       {
1631         spectraName += "-";
1632         spectraName += flavour;
1633       }
1634       spectra = new AliAnalysisMuMuSpectra(spectraName.Data());
1635     }
1636
1637     spectra->AdoptResult(*bin,r);
1638
1639     if ( IsSimulation() )
1640     {
1641       SetNofInputParticles(*r);
1642     }
1643
1644
1645   } // loop on bins
1646
1647   delete fitTypeArray;
1648   delete bins;
1649
1650   return spectra;
1651 }
1652
1653 //_____________________________________________________________________________
1654 std::vector<Double_t>
1655 AliAnalysisMuMu::GetMCCB2Tails(const AliAnalysisMuMuBinning::Range& bin) const
1656 {
1657   /// Get the tails from the associated simulation
1658
1659   std::vector<Double_t> par;
1660
1661   if (!SIM())
1662   {
1663     AliError("Cannot get MC tails without an associated simulation file !");
1664     return par;
1665   }
1666
1667
1668   AliAnalysisMuMuSpectra* s = static_cast<AliAnalysisMuMuSpectra*>(SIM()->GetSpectra(bin.Quantity().Data(),bin.Flavour().Data()));
1669
1670   if (!s)
1671   {
1672     AliError(Form("Could not find spectra %s,%s for associated simulation",bin.Quantity().Data(),bin.Flavour().Data()));
1673     fAssociatedSimulation->MC()->Print("*:Ali*");
1674     return par;
1675   }
1676   else
1677   {
1678     AliDebug(1,Form("AliAnalysisMuMuSpectra* s = reinterpret_cast<AliAnalysisMuMuSpectra*>(%p)",s));
1679   }
1680
1681   AliAnalysisMuMuResult* r = s->GetResultForBin(bin);
1682
1683   if ( r )
1684   {
1685     AliAnalysisMuMuJpsiResult* r1 = dynamic_cast<AliAnalysisMuMuJpsiResult*>(r->SubResult("JPSI:1"));
1686     if  (r1)
1687     {
1688       TF1* func = static_cast<TF1*>(r1->Minv()->GetListOfFunctions()->FindObject("fitTotal"));
1689       if (func)
1690       {
1691         par.push_back(func->GetParameter("alphaLow"));
1692         par.push_back(func->GetParameter("nLow"));
1693         par.push_back(func->GetParameter("alphaUp"));
1694         par.push_back(func->GetParameter("nUp"));
1695       }
1696     }
1697   }
1698
1699   return par;
1700 }
1701
1702 //_____________________________________________________________________________
1703 ULong64_t AliAnalysisMuMu::GetTriggerScalerCount(const char* triggerList, Int_t runNumber)
1704 {
1705   // Get the expected (from OCDB scalers) trigger count
1706
1707   AliAnalysisTriggerScalers ts(runNumber,fgOCDBPath.Data());
1708
1709   TObjArray* triggers = TString(triggerList).Tokenize(",");
1710   TObjString* trigger;
1711   TIter next(triggers);
1712   ULong64_t n(0);
1713
1714   while ( ( trigger = static_cast<TObjString*>(next()) ) )
1715   {
1716     AliAnalysisTriggerScalerItem* item = ts.GetTriggerScaler(runNumber,"L2A",trigger->String().Data());
1717     if (item)
1718     {
1719       n += item->Value();
1720     }
1721     delete item;
1722   }
1723   delete triggers;
1724
1725   return n;
1726 }
1727
1728 //_____________________________________________________________________________
1729 AliAnalysisMuMuSpectra* AliAnalysisMuMu::GetSpectra(const char* what, const char* flavour) const
1730 {
1731   /// get a given spectra
1732
1733   TString swhat(what);
1734   TString sflavour(flavour);
1735   swhat.ToUpper();
1736   sflavour.ToUpper();
1737
1738   TString spectraName(Form("/PSALL/%s/PP/%s/PSI-%s",
1739                            First(fDimuonTriggers).Data(),
1740                            First(fPairSelectionList).Data(),
1741                            swhat.Data()));
1742
1743   if (sflavour.Length()>0)
1744   {
1745     spectraName += "-";
1746     spectraName += sflavour.Data();
1747   }
1748
1749   if (IsSimulation())
1750   {
1751     spectraName.ReplaceAll("PSALL",fgDefaultEventSelectionForSimulations.Data());
1752     spectraName.ReplaceAll(First(fgDefaultDimuonTriggers).Data(),fgDefaultDimuonTriggerForSimulations.Data());
1753   }
1754
1755   return SPECTRA(spectraName.Data());
1756 }
1757
1758 //_____________________________________________________________________________
1759 UInt_t AliAnalysisMuMu::GetSum(AliCounterCollection& cc, const char* triggerList,
1760                                const char* eventSelection, Int_t runNumber)
1761 {
1762   TObjArray* ktrigger = cc.GetKeyWords("trigger").Tokenize(",");
1763   TObjArray* kevent = cc.GetKeyWords("event").Tokenize(",");
1764   TObjArray* a = TString(triggerList).Tokenize(" ");
1765   TIter next(a);
1766   TObjString* str;
1767
1768   UInt_t n(0);
1769
1770   TString sEventSelection(eventSelection);
1771   sEventSelection.ToUpper();
1772
1773   if ( kevent->FindObject(sEventSelection.Data()) )
1774   {
1775     while ( ( str = static_cast<TObjString*>(next()) ) )
1776     {
1777       if ( ktrigger->FindObject(str->String().Data()) )
1778       {
1779         if ( runNumber < 0 )
1780         {
1781           n +=  static_cast<UInt_t>(cc.GetSum(Form("trigger:%s/event:%s",str->String().Data(),eventSelection)));
1782         }
1783         else
1784         {
1785           n +=  static_cast<UInt_t>(cc.GetSum(Form("trigger:%s/event:%s/run:%d",str->String().Data(),eventSelection,runNumber)));
1786         }
1787       }
1788     }
1789   }
1790
1791   delete a;
1792   delete ktrigger;
1793   delete kevent;
1794   return n;
1795 }
1796
1797 //_____________________________________________________________________________
1798 Bool_t
1799 AliAnalysisMuMu::GetCollections(const char* rootfile,
1800                                 AliMergeableCollection*& mc,
1801                                 AliCounterCollection*& cc,
1802                                 AliAnalysisMuMuBinning*& bin,
1803                                 std::set<int>& runnumbers)
1804 {
1805   mc = 0x0;
1806   cc = 0x0;
1807   bin = 0x0;
1808
1809   TFile* f = static_cast<TFile*>(gROOT->GetListOfFiles()->FindObject(rootfile));
1810   if (!f)
1811   {
1812     f = TFile::Open(rootfile);
1813   }
1814
1815   if ( !f || f->IsZombie() )
1816   {
1817     return kFALSE;
1818   }
1819
1820   f->GetObject("MC",mc);
1821   f->GetObject("CC",cc);
1822
1823   TIter next(f->GetListOfKeys());
1824   TKey* key;
1825
1826   while ( ( key = static_cast<TKey*>(next())) && !bin )
1827   {
1828     if ( strcmp(key->GetClassName(),"AliAnalysisMuMuBinning")==0 )
1829     {
1830       bin = dynamic_cast<AliAnalysisMuMuBinning*>(key->ReadObj());
1831     }
1832   }
1833
1834   if (!mc || !cc)
1835   {
1836     AliErrorClass("Old file. Please upgrade it!");
1837
1838     return kFALSE;
1839   }
1840
1841   // get run list
1842   TObjArray* runs = cc->GetKeyWords("run").Tokenize(",");
1843   runs->Sort();
1844   TIter nextRun(runs);
1845   TObjString* srun;
1846
1847   runnumbers.clear();
1848
1849   while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
1850   {
1851     runnumbers.insert(srun->String().Atoi());
1852   }
1853
1854   delete runs;
1855
1856   return kTRUE;
1857 }
1858
1859 //_____________________________________________________________________________
1860 Bool_t AliAnalysisMuMu::IsSimulation() const
1861 {
1862   // whether or not we have MC information
1863
1864   if (!fMergeableCollection) return kFALSE;
1865
1866   return ( fMergeableCollection->Histo(Form("/INPUT/%s/MinvUS",fgDefaultEventSelectionForSimulations.Data())) != 0x0 );
1867 }
1868
1869 //_____________________________________________________________________________
1870 Int_t
1871 AliAnalysisMuMu::Jpsi(const char* what, const char* binningFlavour)
1872 {
1873   // Fit the J/psi (and psiprime) peaks for the triggers in fDimuonTriggers list
1874   // what="integrated" => fit only fully integrated MinvUS
1875   // what="pt" => fit MinvUS in pt bins
1876   // what="y" => fit MinvUS in y bins
1877   // what="pt,y" => fit MinvUS in (pt,y) bins
1878
1879   TStopwatch timer;
1880
1881   if (!fMergeableCollection)
1882   {
1883     AliError("No mergeable collection. Consider Upgrade()");
1884     return 0;
1885   }
1886
1887   Int_t nfits(0);
1888
1889   TObjArray* triggerArray = fDimuonTriggers.Tokenize(",");
1890   TObjArray* eventTypeArray = fEventSelectionList.Tokenize(",");
1891   TObjArray* pairCutArray = fPairSelectionList.Tokenize(",");
1892   TObjArray* whatArray = TString(what).Tokenize(",");
1893
1894   TIter nextTrigger(triggerArray);
1895   TIter nextEventType(eventTypeArray);
1896   TIter nextPairCut(pairCutArray);
1897   TIter nextWhat(whatArray);
1898
1899   TObjString* trigger;
1900   TObjString* eventType;
1901   TObjString* pairCut;
1902   TObjString* swhat;
1903
1904   while ( ( swhat = static_cast<TObjString*>(nextWhat()) ) )
1905   {
1906     AliAnalysisMuMuBinning* binning(0x0);
1907
1908     if ( fBinning && swhat->String().Length() > 0 )
1909     {
1910       binning = fBinning->Project("psi",swhat->String().Data(),binningFlavour);
1911     }
1912     else
1913     {
1914       binning = new AliAnalysisMuMuBinning;
1915       binning->AddBin("psi",swhat->String().Data());
1916     }
1917
1918     StdoutToAliDebug(1,std::cout << "++++++++++++ swhat=" << swhat->String().Data() << std::endl;);
1919
1920     if (!binning)
1921     {
1922       AliError("oups. binning is NULL");
1923       continue;
1924     }
1925
1926     StdoutToAliDebug(1,binning->Print(););
1927
1928     nextTrigger.Reset();
1929
1930     while ( ( trigger = static_cast<TObjString*>(nextTrigger())) )
1931     {
1932       AliDebug(1,Form("TRIGGER %s",trigger->String().Data()));
1933
1934       nextEventType.Reset();
1935
1936       while ( ( eventType = static_cast<TObjString*>(nextEventType())) )
1937       {
1938         AliDebug(1,Form("--EVENTTYPE %s",eventType->String().Data()));
1939
1940         nextPairCut.Reset();
1941
1942         while ( ( pairCut = static_cast<TObjString*>(nextPairCut())) )
1943         {
1944           AliDebug(1,Form("----PAIRCUT %s",pairCut->String().Data()));
1945
1946           AliDebug(1,"----Fitting...");
1947
1948           AliAnalysisMuMuSpectra* spectra = FitParticle("psi",
1949                                                   trigger->String().Data(),
1950                                                   eventType->String().Data(),
1951                                                   pairCut->String().Data(),
1952                                                   "PP",
1953                                                   *binning);
1954
1955           AliDebug(1,Form("----fitting done spectra = %p",spectra));
1956
1957           if ( spectra )
1958           {
1959             ++nfits;
1960
1961             TString id(Form("/%s/%s/PP/%s",eventType->String().Data(),
1962                             trigger->String().Data(),
1963                             pairCut->String().Data()));
1964
1965             TObject* o = fMergeableCollection->GetObject(id.Data(),spectra->GetName());
1966
1967             AliDebug(1,Form("----nfits=%d id=%s o=%p",nfits,id.Data(),o));
1968
1969             if (o)
1970             {
1971               AliWarning(Form("Replacing %s/%s",id.Data(),spectra->GetName()));
1972               fMergeableCollection->Remove(Form("%s/%s",id.Data(),spectra->GetName()));
1973             }
1974
1975             fMergeableCollection->Adopt(id.Data(),spectra);
1976
1977             StdoutToAliDebug(1,spectra->Print(););
1978           }
1979         }
1980       }
1981     }
1982   }
1983
1984   delete whatArray;
1985   delete triggerArray;
1986   delete eventTypeArray;
1987   delete pairCutArray;
1988
1989   StdoutToAliDebug(1,timer.Print(););
1990
1991   if (nfits)
1992   {
1993     Update();
1994 //    ReOpen(fFilename,"UPDATE");
1995 //    fMergeableCollection->Write("MC",TObjArray::kOverwrite);// | TObjArray::kSingleKey);
1996 //    ReOpen(fFilename,"READ");
1997   }
1998
1999
2000   return nfits;
2001
2002 }
2003
2004 //_____________________________________________________________________________
2005 void AliAnalysisMuMu::LatexOutputFnorm(const char* filelist, const char* subresultnames, Bool_t rms)
2006 {
2007   /// Make a LaTeX output of the Fnorm factors for each file in filelist
2008
2009   TObjArray* periods = ReadFileList(filelist);
2010
2011   if (!periods || periods->IsEmpty() ) return;
2012
2013   TIter next(periods);
2014
2015   TObjArray* subresults = TString(subresultnames).Tokenize(",");
2016   TIter nextSub(subresults);
2017
2018   Int_t ic(0);
2019   Int_t iclast = subresults->GetLast();
2020
2021
2022   std::cout << "\\begin{tabular}{l|r|r|r|r}" << std::endl;
2023
2024   std::cout << "Period & ";
2025
2026   TObjString* rname;
2027
2028
2029   while ( ( rname = static_cast<TObjString*>(nextSub())) )
2030   {
2031     std::cout << rname->String().Data();
2032
2033     if ( ic != iclast )
2034     {
2035       std::cout << " & ";
2036     }
2037     ++ic;
2038   }
2039
2040   std::cout << "\\\\" << std::endl << "\\hline" << std::endl;
2041
2042   for ( Int_t i = 0; i <= periods->GetLast(); ++i )
2043   {
2044     TString period("unknown");
2045
2046     TString filename(static_cast<TObjString*>(periods->At(i))->String());
2047
2048     Int_t dummy(0);
2049
2050     if (!DecodeFileName(filename,period,dummy,dummy,dummy))
2051     {
2052       continue;
2053     }
2054
2055     if ( gSystem->AccessPathName(filename) )
2056     {
2057       AliErrorClass(Form("Could not find file %s. Skipping it.",filename.Data()));
2058       continue;
2059
2060     }
2061
2062     AliAnalysisMuMu m(filename.Data());
2063
2064     AliMergeableCollection* norm = m.MC()->Project("/FNORM/RESULTS/");
2065
2066     AliAnalysisMuMuResult* r = static_cast<AliAnalysisMuMuResult*>(norm->GetObject("Fnorm")->Clone());
2067
2068     nextSub.Reset();
2069
2070     if (!r)
2071     {
2072       AliErrorClass(Form("Could not get result for file %s.",filename.Data()));
2073       continue;
2074     }
2075
2076     std::cout << period << " & ";
2077
2078     next.Reset();
2079
2080     ic = 0;
2081
2082     while ( ( rname = static_cast<TObjString*>(nextSub())) )
2083     {
2084       TString name(rname->String());
2085       TString var("Fnorm");
2086
2087       if ( name.BeginsWith("RelDif"))
2088       {
2089         r = static_cast<AliAnalysisMuMuResult*>(norm->GetObject("RelDif"));
2090
2091       }
2092       else if ( name.BeginsWith("Correction"))
2093       {
2094         r = static_cast<AliAnalysisMuMuResult*>(norm->GetObject("Correction"));
2095       }
2096       else
2097       {
2098         r = static_cast<AliAnalysisMuMuResult*>(norm->GetObject("Fnorm"));
2099       }
2100       r->Exclude("*");
2101       r->Include(name.Data());
2102
2103       if ( rms )
2104       {
2105         std::cout << Form(" $%7.2f (%5.2f \%%)$ ",
2106                           r->GetRMS(var.Data()),
2107                           100.0*r->GetRMS(var.Data())/r->GetValue(var.Data()));
2108
2109       }
2110       else
2111       {
2112         std::cout << Form(" $%7.2f \\pm %5.3f$  ",r->GetValue(var.Data()),
2113                           r->GetErrorStat(var.Data()));
2114       }
2115
2116       if ( ic != iclast )
2117       {
2118         std::cout << "&";
2119       }
2120
2121       ++ic;
2122     }
2123
2124     if ( i != periods->GetLast() )
2125     {
2126       std::cout << "\\\\";
2127     }
2128
2129     std::cout << std::endl;
2130
2131     delete norm;
2132   }
2133
2134
2135   std::cout << "\\end{tabular}" << std::endl;
2136
2137
2138
2139   delete periods;
2140 }
2141
2142 //_____________________________________________________________________________
2143 void AliAnalysisMuMu::PlotBackgroundEvolution(const char* gfile, const char* triggerList, Double_t ymax, Bool_t fillBoundaries)
2144 {
2145   // plot the graphs found in the file (and generated using the ComputeBackgroundEvolution() method)
2146
2147   TFile* f = TFile::Open(ExpandPathName(gfile).Data());
2148
2149   if ( !f || !f->IsOpen() )
2150   {
2151     return;
2152   }
2153
2154   SetColorScheme();
2155
2156
2157   TCanvas* c = new TCanvas("background-evolution","background-evolution");
2158
2159   c->Draw();
2160
2161   TLegend* l = new TLegend(0.4,0.6,0.97,0.97);
2162   l->SetFillColor(0);
2163   l->SetTextColor(AliAnalysisMuMu::kBlue);
2164   l->SetLineColor(AliAnalysisMuMu::kBlue);
2165
2166   TObjArray* triggers = TString(triggerList).Tokenize(",");
2167
2168   gStyle->SetOptTitle(0);
2169
2170   TObjString* str(0x0);
2171   TIter next(triggers);
2172   Int_t i(0);
2173   Int_t run1(99999999);
2174   Int_t run2(0);
2175
2176   std::set<int> runnumbers;
2177
2178   while ( ( str = static_cast<TObjString*>(next()) ) )
2179   {
2180     TGraph* g = static_cast<TGraph*>(f->Get(Form("mu_%s",str->String().Data())));
2181     if (!g) continue;
2182     for ( Int_t ir = 0; ir < g->GetN(); ++ir )
2183     {
2184       Int_t run = TMath::Nint(g->GetX()[ir]);
2185       runnumbers.insert(run);
2186       run1 = TMath::Min(run1,run);
2187       run2 = TMath::Max(run2,run);
2188     }
2189   }
2190
2191   AliInfoClass(Form("run1 %d run2 %d",run1,run2));
2192
2193   Double_t ymin(0);
2194
2195   TH2* hframe = new TH2F("hframe","hframe",run2-run1+1,run1,run2,100,ymin,ymax);
2196   hframe->Draw();
2197   hframe->GetXaxis()->SetNoExponent();
2198   hframe->GetYaxis()->SetTitle("Background percentage");
2199
2200   if (fillBoundaries)
2201   {
2202     AliAnalysisTriggerScalers ts(runnumbers,fgOCDBPath.Data());
2203     ts.DrawFills(ymin,ymax);
2204   }
2205
2206   next.Reset();
2207
2208   while ( ( str = static_cast<TObjString*>(next()) ) )
2209   {
2210     TGraph* g = static_cast<TGraph*>(f->Get(Form("mu_%s",str->String().Data())));
2211     if (!g)
2212     {
2213       AliErrorClass(Form("Graph mu_%s not found",str->String().Data()));
2214       continue;
2215     }
2216
2217     Int_t color(i+1);
2218
2219     if (i==0) color = AliAnalysisMuMu::kBlue;
2220     if (i==1) color = AliAnalysisMuMu::kOrange;
2221
2222     g->SetLineColor(color);
2223     g->SetMarkerColor(color);
2224     g->SetMarkerStyle(20+i);
2225
2226     g->Draw("LP");
2227
2228     TLegendEntry* le = l->AddEntry(g,str->String().Data(),"lp");
2229     le->SetTextColor(color);
2230
2231     g->GetYaxis()->SetTitleColor(AliAnalysisMuMu::kBlue);
2232     g->GetXaxis()->SetTitleColor(AliAnalysisMuMu::kBlue);
2233     //    g->Print();
2234
2235     ++i;
2236   }
2237
2238   hframe->Draw("sameaxis");
2239
2240   l->Draw();
2241 }
2242
2243 //_____________________________________________________________________________
2244 void
2245 AliAnalysisMuMu::PlotJpsiEvolution(const char* resultFile, const char* triggerList, Bool_t fillBoundaries,
2246                                    const char* efficiencyFile, Bool_t simulation)
2247 {
2248   /// Will plot the Jpsi rate (or AccxEff if simulation=kTRUE) evolution.
2249   /// (JpsiRate is the number of Jpsi divided by the number of triggers)
2250
2251   std::map<int, float> efficiencies;
2252
2253   if ( efficiencyFile && strlen(efficiencyFile) > 0 )
2254   {
2255     std::ifstream in(gSystem->ExpandPathName(efficiencyFile));
2256     if (!in.bad())
2257     {
2258       char line[1024];
2259       int run;
2260       float eff;
2261       float dummy,errorl,errorh;
2262       int idummy;
2263       while ( in.getline(line,1023,'\n') )
2264       {
2265         sscanf(line,"%d, x[%d]=%f, y[%d]=%f, exl[%d]=%f, exh[%d]=%f, eyl[%d]=%f, eyh[%d]=%f",
2266                &run,&idummy,&dummy,&idummy,&eff,&idummy,&dummy,&idummy,&dummy,&idummy,&errorl,&idummy,&errorh);
2267
2268         AliDebugClass(1,Form("%09d %8.6f +- %8.6f ",run,eff,errorl+errorh));
2269
2270         efficiencies[run] = eff;
2271       }
2272     }
2273   }
2274
2275   TFile* f = TFile::Open(gSystem->ExpandPathName(resultFile));
2276
2277   std::set<int> runnumbers;
2278
2279   TMap* m = static_cast<TMap*>(f->Get("rbr"));
2280
2281   TIter next(m);
2282   TObjString* str;
2283
2284   TObjArray files;
2285   files.SetOwner(kTRUE);
2286
2287   while ( ( str = static_cast<TObjString*>(next())) )
2288   {
2289     files.Add(new TObjString(str->String()));
2290   }
2291
2292   files.Sort();
2293
2294   std::map<std::string, std::vector<float> > x_jpsirate;
2295   std::map<std::string, std::vector<float> > y_jpsirate;
2296   std::map<std::string, std::vector<float> > xerr_jpsirate;
2297   std::map<std::string, std::vector<float> > yerr_jpsirate;
2298
2299   TIter nextTrigger(TString(triggerList).Tokenize(","));
2300   TObjString* trigger(0x0);
2301
2302   int runMin(100000000);
2303   int runMax(0);
2304
2305   TIter nextFile(&files);
2306
2307   Double_t ymin(TMath::Limits<double>::Max());
2308   Double_t ymax(TMath::Limits<double>::Min());
2309
2310
2311   while ( ( trigger = static_cast<TObjString*>(nextTrigger())))
2312   {
2313     Double_t sumw(0);
2314     Double_t n(0);
2315
2316     TString triggerClass(trigger->String());
2317
2318     nextFile.Reset();
2319
2320     while ( ( str = static_cast<TObjString*>(nextFile())) )
2321     {
2322       TObjArray* a = static_cast<TObjArray*>(m->GetValue(str->String().Data()));
2323       if (!a) continue;
2324       AliAnalysisMuMuJpsiResult* r = static_cast<AliAnalysisMuMuJpsiResult*>(a->FindObject(triggerClass.Data()));
2325       if (!r) continue;
2326
2327       TString period;
2328       int aodtrain,esdpass,runnumber;
2329
2330       if ( DecodeFileName(str->String().Data(),period,esdpass,aodtrain,runnumber) )
2331       {
2332         runnumbers.insert(runnumber);
2333
2334         runMin = TMath::Min(runMin,runnumber);
2335         runMax = TMath::Max(runMax,runnumber);
2336
2337         x_jpsirate[triggerClass.Data()].push_back(runnumber);
2338         xerr_jpsirate[triggerClass.Data()].push_back(0.5);
2339
2340         Double_t y(0.0);
2341         Double_t yerr(0.0);
2342         TString what("RateJpsi");
2343         if ( simulation )
2344         {
2345           what = "AccEffJpsi";
2346         }
2347
2348         if ( TMath::Finite(r->GetValue("SigmaJpsi")) && r->NofTriggers() > 10 )
2349         {
2350           y = 100*r->GetValue(what.Data());
2351           yerr = 100*r->GetErrorStat(what.Data());
2352
2353           if  (!efficiencies.empty() )
2354           {
2355             if (efficiencies.count(runnumber))
2356             {
2357               y /= ( efficiencies[runnumber] );
2358             }
2359             else
2360             {
2361               continue;
2362             }
2363           }
2364
2365           ymin = TMath::Min(ymin,y);
2366           ymax = TMath::Max(ymax,y);
2367
2368           sumw += y*r->NofTriggers();
2369           n += r->NofTriggers();
2370         }
2371
2372         y_jpsirate[triggerClass.Data()].push_back(y);
2373         yerr_jpsirate[triggerClass.Data()].push_back(yerr);
2374       }
2375     }
2376
2377     AliInfoClass(Form("Trigger %30s ponderated mean is %7.2f",trigger->String().Data(),sumw/n));
2378   }
2379
2380   delete f;
2381
2382   TString canvasName("cJpsiRateEvolution");
2383
2384   if ( !efficiencies.empty() )
2385   {
2386     canvasName += "Corr";
2387
2388   }
2389   TCanvas* c = new TCanvas(canvasName.Data(),canvasName.Data());
2390
2391   c->Draw();
2392
2393   Int_t nbins = runnumbers.size();
2394   Int_t xmin(runMin);
2395   Int_t xmax(runMax);
2396
2397   if ( CompactGraphs() )
2398   {
2399     xmin = 0;
2400     xmax = nbins-1;
2401   }
2402
2403   TH2* h = new TH2F("h",Form("h;RunNumber;%s",(simulation ? "AccxEff (%)":"J/#psi per CMUL (%)")),
2404                     nbins,xmin,xmax,100,ymin,ymax*1.2);
2405
2406   gStyle->SetOptTitle(0);
2407   gStyle->SetOptStat(0);
2408
2409   if (!CompactGraphs())
2410   {
2411     h->GetXaxis()->SetNoExponent();
2412   }
2413   else
2414   {
2415     std::set<int>::const_iterator it;
2416     int i(0);
2417
2418     for ( it = runnumbers.begin(); it != runnumbers.end(); ++it )
2419     {
2420       h->GetXaxis()->SetBinLabel(i,Form("%d",*it));
2421       ++i;
2422     }
2423     h->GetXaxis()->SetNdivisions(1,kFALSE);
2424
2425   }
2426
2427   h->Draw();
2428
2429   if (fillBoundaries)
2430   {
2431     AliAnalysisTriggerScalers ts(runnumbers,fgOCDBPath);
2432     ts.DrawFills(ymin,ymax);
2433   }
2434
2435   h->Draw("sameaxis");
2436
2437   //c->RedrawAxis("g");
2438
2439   nextTrigger.Reset();
2440
2441   int i(0);
2442   int color[] = { 2,1,4,5,6 };
2443   int marker[] = { 20,23,25,21,22 };
2444
2445   while ( ( trigger = static_cast<TObjString*>(nextTrigger())))
2446   {
2447     std::vector<float>& x = x_jpsirate[trigger->String().Data()];
2448     std::vector<float>& y = y_jpsirate[trigger->String().Data()];
2449     std::vector<float>& xerr = xerr_jpsirate[trigger->String().Data()];
2450     std::vector<float>& yerr = yerr_jpsirate[trigger->String().Data()];
2451
2452     TGraphErrors* g = new TGraphErrors(x.size(),&x[0],&y[0],&xerr[0],&yerr[0]);
2453
2454     g->SetLineColor(1);//color[i]);
2455     g->SetMarkerColor(color[i]);
2456     g->SetMarkerStyle(marker[i]);
2457     g->SetMarkerSize(0.7);
2458     g->GetXaxis()->SetNoExponent();
2459
2460     if ( CompactGraphs() )
2461     {
2462       AliAnalysisMuMuGraphUtil::Compact(*g);
2463     }
2464
2465     g->Draw("P");
2466     TString gname(trigger->String());
2467     gname.ReplaceAll("-","_");
2468     g->SetName(gname.Data());
2469 //    g->Print();
2470
2471     Double_t m2 = g->GetMean(2);
2472
2473     TLine* line = new TLine(runMin,m2,runMax,m2);
2474     line->SetLineColor(color[i]);
2475     line->Draw();
2476
2477     AliInfoClass(Form("TRIGGER %s MEAN %7.2f",trigger->String().Data(),m2));
2478     ++i;
2479   }
2480
2481
2482 }
2483
2484 //_____________________________________________________________________________
2485 TGraph* AliAnalysisMuMu::PlotEventSelectionEvolution(const char* trigger1, const char* event1,
2486                                                      const char* trigger2, const char* event2,
2487                                                      Bool_t drawFills,
2488                                                      Bool_t asRejection) const
2489 {
2490   if (!CC()) return 0x0;
2491
2492   const std::set<int>& runnumbers = RunNumbers();
2493
2494   TGraphErrors* g = new TGraphErrors(runnumbers.size());
2495
2496   std::set<int>::const_iterator it;
2497   Int_t i(0);
2498
2499   Double_t ymin(TMath::Limits<double>::Max());
2500   Double_t ymax(TMath::Limits<double>::Min());
2501
2502   for ( it = runnumbers.begin(); it != runnumbers.end(); ++it )
2503   {
2504     Int_t runNumber = *it;
2505     Double_t n = CC()->GetSum(Form("trigger:%s/event:%s/run:%d",trigger1,event1,runNumber));
2506     Double_t d = CC()->GetSum(Form("trigger:%s/event:%s/run:%d",trigger2,event2,runNumber));
2507     if (n>0 && d>0)
2508     {
2509       Double_t y = n/d;
2510
2511       if ( fCorrectionPerRun )
2512       {
2513         Double_t xcorr,ycorr;
2514         fCorrectionPerRun->GetPoint(i,xcorr,ycorr); // note that the fact that xcorr==runNumber has been checked by the SetCorrectionPerRun method
2515         y *= ycorr;
2516         // FIXME: should get the correction error here
2517       }
2518
2519       if ( asRejection ) y = 100*(1.0 - y);
2520       ymin = TMath::Min(ymin,y);
2521       ymax = TMath::Max(ymax,y);
2522       Double_t yerr = y*AliAnalysisMuMuResult::ErrorAB(n,TMath::Sqrt(n),d,TMath::Sqrt(d));
2523       g->SetPoint(i,runNumber,y);
2524       g->SetPointError(i,0.5,yerr);
2525
2526       ++i;
2527     }
2528
2529   }
2530
2531   TH2* hframe = new TH2F(Form("%s %s-%s / %s-%s",(asRejection ? "1 - ":""),trigger1,event1,trigger2,event2),
2532                          Form("%s %s-%s / %s-%s",(asRejection ? "1 - ":""),trigger1,event1,trigger2,event2),
2533                          runnumbers.size()+50,
2534                          *(runnumbers.begin())-25,
2535                          *(runnumbers.rbegin())+25,100,0,ymax*1.3);
2536
2537   gStyle->SetOptStat(0);
2538
2539   hframe->Draw();
2540
2541   hframe->GetXaxis()->SetNoExponent();
2542
2543   hframe->GetYaxis()->SetTitle(asRejection ? "Rejection (%)" : "Ratio");
2544
2545   g->Set(i);
2546   g->SetTitle(Form("%s %s-%s / %s-%s",(asRejection ? "1 - ":""),trigger1,event1,trigger2,event2));
2547   g->GetXaxis()->SetNoExponent();
2548   g->Draw("lp");
2549
2550   AliAnalysisTriggerScalers ts(RunNumbers(),fgOCDBPath.Data());
2551
2552   if ( drawFills )
2553   {
2554     ts.DrawFills(ymin,ymax);
2555     g->Draw("lp");
2556   }
2557
2558
2559   std::map<std::string, std::pair<int,int> > periods;
2560
2561   ts.GetLHCPeriodBoundaries(periods);
2562
2563   TLegend* legend = new TLegend(0.15,0.82,0.90,0.92);
2564   legend->SetFillColor(0);
2565   Int_t n(0);
2566
2567
2568   for ( std::map<std::string, std::pair<int,int> >::const_iterator pit = periods.begin(); pit != periods.end(); ++pit )
2569   {
2570     std::string period = pit->first;
2571     int run1 = (pit->second).first;
2572     int run2 = (pit->second).second;
2573     int nruns(0);
2574     for ( std::set<int>::const_iterator rit = RunNumbers().begin(); rit != RunNumbers().end(); ++ rit )
2575     {
2576       if ( (*rit) >= run1 && (*rit) <= run2 )
2577       {
2578         ++nruns;
2579       }
2580     }
2581     AliInfo(Form("Period %s runs %6d-%6d ; %d actual runs",period.c_str(),run1,run2,nruns));
2582
2583     g->Fit("pol0","+Q","",run1,run2);
2584     TF1* func = static_cast<TF1*>(g->GetListOfFunctions()->Last());
2585     if (func)
2586     {
2587       func->SetLineColor(2+n);
2588       legend->AddEntry(func,Form("%s %5.2f #pm %5.2f %s (rel. error %5.2f %%)",period.c_str(),func->GetParameter(0),func->GetParError(0),
2589                                  (asRejection ? "%":""),100*func->GetParError(0)/func->GetParameter(0)));
2590       ++n;
2591     }
2592   }
2593
2594   legend->SetNColumns(3);
2595
2596   Double_t mean = TMath::Mean(g->GetN(),g->GetY());
2597   Double_t rms = TMath::RMS(g->GetN(),g->GetY());
2598
2599   legend->AddEntry("",Form("Mean %5.2f RMS %5.2f (%5.2f %%)",mean,rms,(mean) ? 100.0*rms/mean : 0.0),"");
2600
2601   legend->Draw();
2602
2603   return g;
2604 }
2605
2606 //_____________________________________________________________________________
2607 void AliAnalysisMuMu::Print(Option_t* opt) const
2608 {
2609     /// printout
2610   std::cout << "Reading from file : " << fFilename.Data() << std::endl;
2611   std::cout << "List of dimuon triggers to consider : " << DimuonTriggerList() << std::endl;
2612   std::cout << "List of   muon triggers to consider : " << MuonTriggerList() << std::endl;
2613   std::cout << "List of     MB triggers to consider : " << MinbiasTriggerList() << std::endl;
2614   std::cout << "Event selection list : " << EventSelectionList() << std::endl;
2615   std::cout << "Pair  selection list : " << PairSelectionList() << std::endl;
2616
2617   std::cout << RunNumbers().size() << " runs";
2618   if ( fCorrectionPerRun )
2619   {
2620     std::cout << " with correction factors";
2621   }
2622   std::cout << std::endl;
2623   Int_t i(0);
2624   for ( std::set<int>::const_iterator it = RunNumbers().begin(); it != RunNumbers().end(); ++it )
2625   {
2626     std::cout << (*it);
2627     if ( fCorrectionPerRun )
2628     {
2629       std::cout << Form("(%e)",fCorrectionPerRun->GetY()[i]);
2630     }
2631     std::cout << ",";
2632     ++i;
2633   }
2634   std::cout << std::endl;
2635
2636   TString sopt(opt);
2637   sopt.ToUpper();
2638
2639   if ( sopt.Contains("BIN") && BIN() )
2640   {
2641     std::cout << "Binning : " << std::endl;
2642     TString topt(sopt);
2643     topt.ReplaceAll("BIN","");
2644     BIN()->Print(topt.Data());
2645   }
2646   if ( sopt.Contains("MC") && MC() )
2647   {
2648     TString topt(sopt);
2649     topt.ReplaceAll("MC","");
2650     MC()->Print(topt.Data());
2651   }
2652   if ( sopt.Contains("CC") && CC() )
2653   {
2654     CC()->Print("trigger/event");
2655   }
2656
2657   if ( sopt.Contains("SIZE") )
2658   {
2659     TFile* f = ReOpen(fFilename,"READ");
2660     TIter next(f->GetListOfKeys());
2661     TKey* key;
2662
2663     while ( ( key = static_cast<TKey*>(next()) ) )
2664     {
2665       std::cout << key->GetName() << " " << key->GetNbytes() << " " << key->GetObjlen() << std::endl;
2666     }
2667   }
2668 }
2669
2670 //_____________________________________________________________________________
2671 TObjArray*
2672 AliAnalysisMuMu::ReadFileList(const char* filelist)
2673 {
2674   //
2675   // read the filelist and try to order it by runnumber (or periods)
2676   //
2677   // filelist can either be a real filelist (i.e. a text file containing
2678   // root filenames) or a root file itself.
2679   //
2680
2681   char line[1024];
2682
2683   TObjArray* files = new TObjArray;
2684   files->SetOwner(kTRUE);
2685
2686   TString sfilelist(ExpandPathName(filelist));
2687
2688   if ( sfilelist.EndsWith(".root") )
2689   {
2690     files->Add(new TObjString(sfilelist.Data()));
2691     return files;
2692   }
2693
2694   std::set<std::string> sorting;
2695   std::map<std::string,std::string> filemap;
2696
2697   std::ifstream in(sfilelist.Data());
2698
2699   TString period;
2700   int aodtrain,esdpass,runnumber;
2701
2702   while ( in.getline(line,1022,'\n') )
2703   {
2704     TString sline(line);
2705     if (sline.BeginsWith("#")) continue;
2706
2707     DecodeFileName(line,period,esdpass,aodtrain,runnumber);
2708
2709     AliDebugClass(1,Form("line %s => period %s esdpass %d aodtrain %d runnumber %09d",
2710                          line,period.Data(),esdpass,aodtrain,runnumber));
2711
2712     TString key(Form("%d",runnumber));
2713
2714     if ( runnumber <= 0 )
2715     {
2716       key = period;
2717     }
2718     sorting.insert(key.Data());
2719     filemap.insert(std::make_pair<std::string,std::string>(key.Data(),line));
2720   }
2721
2722   in.close();
2723
2724   std::set<std::string>::const_iterator it;
2725
2726   for ( it = sorting.begin(); it != sorting.end(); ++it )
2727   {
2728     files->Add(new TObjString(filemap[*it].c_str()));
2729   }
2730   return files;
2731 }
2732
2733 //_____________________________________________________________________________
2734 TFile* AliAnalysisMuMu::ReOpen(const char* filename, const char* mode) const
2735 {
2736   /// Tries to reopen the file with a new mode
2737
2738   TFile* f = static_cast<TFile*>(gROOT->GetListOfFiles()->FindObject(filename));
2739
2740   if (f)
2741   {
2742     delete f;
2743   }
2744
2745   f = TFile::Open(filename,mode);
2746
2747   if ( !f || !f->IsOpen() )
2748   {
2749     AliError(Form("Cannot open file %s in mode %s",filename,mode));
2750     return 0x0;
2751   }
2752
2753   return f;
2754 }
2755
2756 //_____________________________________________________________________________
2757 Int_t
2758 AliAnalysisMuMu::RunNumberFromFileName(const char* filename)
2759 {
2760   TString period;
2761   int esdpass,aodtrain,runnumber;
2762   Bool_t ok = DecodeFileName(filename,period,esdpass,aodtrain,runnumber);
2763   if ( ok ) return runnumber;
2764   return -1;
2765 }
2766
2767 //_____________________________________________________________________________
2768 void AliAnalysisMuMu::SetColorScheme()
2769 {
2770   new TColor(AliAnalysisMuMu::kBlue,4/255.0,44/255.0,87/255.0,"my blue");
2771   new TColor(AliAnalysisMuMu::kOrange,255/255.0,83/255.0,8/255.0,"my orange");
2772   new TColor(AliAnalysisMuMu::kGreen,152/255.0,202/255.0,52/255.0,"my green");
2773
2774   gStyle->SetGridColor(AliAnalysisMuMu::kBlue);
2775
2776   gStyle->SetFrameLineColor(AliAnalysisMuMu::kBlue);
2777   gStyle->SetAxisColor(AliAnalysisMuMu::kBlue,"xyz");
2778   gStyle->SetLabelColor(AliAnalysisMuMu::kBlue,"xyz");
2779
2780   gStyle->SetTitleColor(AliAnalysisMuMu::kBlue);
2781   gStyle->SetTitleTextColor(AliAnalysisMuMu::kBlue);
2782   gStyle->SetLabelColor(AliAnalysisMuMu::kBlue);
2783   gStyle->SetStatTextColor(AliAnalysisMuMu::kBlue);
2784
2785   gStyle->SetOptStat(0);
2786 }
2787
2788 //_____________________________________________________________________________
2789 Bool_t AliAnalysisMuMu::SetCorrectionPerRun(const TGraph& corr, const char* formula)
2790 {
2791     /// Sets the graph used to correct values per run
2792   delete fCorrectionPerRun;
2793   fCorrectionPerRun=0x0;
2794
2795   // check that corr has the same runs as we do
2796
2797   Int_t i(0);
2798
2799   for ( std::set<int>::const_iterator it = RunNumbers().begin(); it != RunNumbers().end(); ++it )
2800   {
2801     Int_t corrRun = TMath::Nint(corr.GetX()[i]);
2802     if (corrRun != *it)
2803     {
2804       AliError(Form("%d-th run mistmatch %d vs %d",i,corrRun,*it));
2805
2806       return kFALSE;
2807     }
2808     ++i;
2809   }
2810
2811   fCorrectionPerRun = new TGraphErrors(corr.GetN());
2812
2813   TFormula* tformula(0x0);
2814   if ( strlen(formula) > 0 )
2815   {
2816     tformula = new TFormula("SetCorrectionPerRunFormula",formula);
2817   }
2818
2819   i = 0;
2820
2821   for ( std::set<int>::const_iterator it = RunNumbers().begin(); it != RunNumbers().end(); ++it )
2822   {
2823     Double_t y = corr.GetY()[i];
2824
2825     if ( tformula )
2826     {
2827       y = tformula->Eval(y);
2828     }
2829     fCorrectionPerRun->SetPoint(i,corr.GetX()[i],y);
2830     ++i;
2831   }
2832
2833   delete formula;
2834
2835   return kTRUE;
2836 }
2837
2838 //_____________________________________________________________________________
2839 void AliAnalysisMuMu::SetNofInputParticles(AliAnalysisMuMuJpsiResult& r)
2840 {
2841   /// Set the "NofInput" variable(s) of one result
2842
2843   TString hname(Form("MinvUS%s",r.Bin().AsString().Data()));
2844
2845   TH1* hinput = fMergeableCollection->Histo("/INPUT/INYRANGE",hname.Data());
2846
2847   if (!hinput)
2848   {
2849     AliError(Form("Got a simulation file where I did not find histogram /INPUT/INYRANGE/%s",hname.Data()));
2850
2851   }
2852   else
2853   {
2854     r.SetNofInputParticles(*hinput);
2855   }
2856 }
2857
2858 //_____________________________________________________________________________
2859 AliAnalysisMuMuSpectra* AliAnalysisMuMu::SPECTRA(const char* fullpath) const
2860 {
2861   /// Shortcut method to get to a spectra
2862   if (!MC()) return 0x0;
2863
2864   return static_cast<AliAnalysisMuMuSpectra*>(MC()->GetObject(fullpath));
2865 }
2866
2867 //_____________________________________________________________________________
2868 void AliAnalysisMuMu::TriggerCountCoverage(const char* triggerList,
2869                                            Bool_t compact,
2870                                            Bool_t orderByTriggerCount)
2871 {
2872   // Give the fraction of triggers (in triggerList) relative
2873   // to what is expected in the scalers
2874
2875   TGrid::Connect("alien://"); // to insure the "Trying to connect to server... message does not pollute our output later on...
2876
2877   AliLog::EType_t oldLevel = static_cast<AliLog::EType_t>(AliLog::GetGlobalLogLevel());
2878
2879   AliLog::SetGlobalLogLevel(AliLog::kFatal);
2880
2881   if (!fMergeableCollection || !fCounterCollection) return;
2882
2883   TObjArray* runs = fCounterCollection->GetKeyWords("run").Tokenize(",");
2884   TIter nextRun(runs);
2885
2886   TObjArray* triggers = fCounterCollection->GetKeyWords("trigger").Tokenize(",");
2887   TIter nextTrigger(triggers);
2888
2889   TObjString* srun;
2890   TObjString* strigger;
2891
2892   TString striggerList(triggerList);
2893
2894   ULong64_t total(0);
2895   ULong64_t totalExpected(0);
2896   TString msg;
2897   std::multimap<ULong64_t,std::string> messages;
2898
2899   while ( ( srun = static_cast<TObjString*>(nextRun()) ) )
2900   {
2901     msg.Form("RUN %09d ",srun->String().Atoi());
2902
2903     if (!compact)
2904     {
2905         msg += "\n";
2906     }
2907
2908     ULong64_t nmax(0);
2909
2910     nextTrigger.Reset();
2911
2912     while ( ( strigger = static_cast<TObjString*>(nextTrigger()) ) )
2913     {
2914       if ( !striggerList.Contains(strigger->String().Data()) )
2915       {
2916         continue;
2917       }
2918       ULong64_t n = TMath::Nint(fCounterCollection->GetSum(Form("trigger:%s/event:%s/run:%d",
2919                                                                 strigger->String().Data(),"ALL",srun->String().Atoi())));
2920
2921       ULong64_t expected = GetTriggerScalerCount(strigger->String().Data(),srun->String().Atoi());
2922
2923
2924       nmax = TMath::Max(n,nmax);
2925
2926       total += n;
2927       totalExpected += expected;
2928
2929       msg += TString::Format("%30s %9lld expected %9lld [%s] ",strigger->String().Data(),n,expected,
2930                              (n>expected ? "!" : " "));
2931
2932       if ( expected > 0 ) {
2933         msg += TString::Format("fraction %5.1f %%",n*100.0/expected);
2934       }
2935
2936       if (!compact)
2937       {
2938         msg += "\n";
2939       }
2940     }
2941     if (nmax>0)
2942     {
2943       if (!orderByTriggerCount)
2944       {
2945         std::cout << msg.Data() << std::endl;
2946       }
2947       else
2948       {
2949         messages.insert(std::make_pair(nmax,static_cast<std::string>(msg.Data())));
2950       }
2951     }
2952   }
2953
2954   std::multimap<ULong64_t,std::string>::const_reverse_iterator it;
2955
2956   ULong64_t current(0);
2957   Int_t n(0);
2958
2959   for ( it = messages.rbegin(); it != messages.rend(); ++it )
2960   {
2961     ++n;
2962     current += it->first;
2963     Double_t percent = 0.0;
2964     if ( total > 0.0 )
2965     {
2966       percent = current*100.0/total;
2967     }
2968     std::cout << Form("%10lld",it->first) << " " << it->second << " percentage of total = " << Form("%7.2f %% %3d",percent,n ) << std::endl;
2969   }
2970
2971   std::cout << Form("--- TOTAL %lld expected %lld fraction %5.1f %%",
2972                     total,totalExpected,totalExpected ? total*100.0/totalExpected : 0.0) << std::endl;
2973
2974
2975   AliLog::SetGlobalLogLevel(oldLevel);
2976   delete triggers;
2977   delete runs;
2978 }
2979
2980 //_____________________________________________________________________________
2981 void AliAnalysisMuMu::UnsetCorrectionPerRun()
2982 {
2983     // drop the correction factors
2984   delete fCorrectionPerRun;
2985   fCorrectionPerRun=0x0;
2986 }
2987
2988 //_____________________________________________________________________________
2989 void AliAnalysisMuMu::Update()
2990 {
2991   /// update the current file with memory
2992
2993   if (!CC() || !MC()) return;
2994
2995   ReOpen(fFilename,"UPDATE");
2996
2997   if (MC())
2998   {
2999     MC()->Write("MC",TObject::kSingleKey|TObject::kOverwrite);
3000   }
3001
3002   ReOpen(fFilename,"READ");
3003
3004   GetCollections(fFilename,fMergeableCollection,fCounterCollection,fBinning,fRunNumbers);
3005 }
3006
3007 //_____________________________________________________________________________
3008 Bool_t AliAnalysisMuMu::Upgrade(const char* filename)
3009 {
3010   /// Upgrade a file
3011   AliAnalysisMuMu m(filename);
3012
3013   return m.Upgrade();
3014 }
3015
3016 //_____________________________________________________________________________
3017 Bool_t AliAnalysisMuMu::Upgrade()
3018 {
3019   /// Upgrade the current file
3020   /// - from single list to one key per object, if needed
3021   /// - from histogramCollection to mergeableCollection, if needed
3022
3023   TFile* f = ReOpen(fFilename,"UPDATE");
3024
3025   TList* list = static_cast<TList*>(f->Get("chist"));
3026
3027   if (list)
3028   {
3029     // really old file where everything was in a single list
3030
3031     AliHistogramCollection* hc = static_cast<AliHistogramCollection*>(list->At(0));
3032     AliCounterCollection* cc = static_cast<AliCounterCollection*>(list->At(1));
3033
3034     AliMergeableCollection* mc = hc->Convert();
3035
3036     f->cd();
3037
3038     mc->Write("MC",TObject::kSingleKey);
3039     cc->Write("CC",TObject::kSingleKey);
3040
3041     f->Delete("chist;*");
3042
3043     f->Write();
3044
3045   }
3046   else
3047   {
3048     AliHistogramCollection* hc = static_cast<AliHistogramCollection*>(f->Get("HC"));
3049
3050     if ( hc )
3051     {
3052       // old file with histogram collection instead of mergeable collection
3053
3054       AliMergeableCollection* mc = hc->Convert();
3055
3056       f->cd();
3057
3058       mc->Write("MC",TObject::kSingleKey);
3059
3060       f->Delete("HC;*");
3061
3062       f->Write();
3063     }
3064   }
3065
3066   delete f;
3067
3068   return kTRUE;
3069 }
3070
3071 //_____________________________________________________________________________
3072 void AliAnalysisMuMu::ComputeFnorm()
3073 {
3074   /// Compute the CMUL to CINT ratio(s)
3075
3076   if (!CC()) return;
3077
3078   MC()->Prune("/FNORM");
3079
3080   AliAnalysisMuMuFnorm computer(*(CC()),AliAnalysisMuMuFnorm::kMUL,fgOCDBPath.Data(),fgIsCompactGraphs);
3081
3082   computer.ComputeFnorm();
3083
3084   AliMergeableCollection* fnorm = computer.DetachMC();
3085
3086   MC()->Attach(fnorm,"/FNORM/");
3087
3088   Update();
3089 }
3090
3091 //_____________________________________________________________________________
3092 AliAnalysisMuMuSpectra* AliAnalysisMuMu::CorrectSpectra(const char* type, const char* flavour)
3093 {
3094   /// Correct one spectra
3095
3096   if (!SIM())
3097   {
3098     AliError("Cannot compute corrected yield without associated MC file !");
3099     return 0x0;
3100   }
3101
3102   const char* accEffSubResultName="COUNTJPSI:1";
3103
3104   AliAnalysisMuMuSpectra* realSpectra = GetSpectra(type,flavour);
3105   AliAnalysisMuMuSpectra* simSpectra = SIM()->GetSpectra(type,flavour);
3106
3107   if ( !realSpectra )
3108   {
3109     AliError("could not get real spectra");
3110     return 0x0;
3111   }
3112
3113   if ( !simSpectra)
3114   {
3115     AliError("could not get sim spectra");
3116     return 0x0;
3117   }
3118
3119   realSpectra->Correct(*simSpectra,"Jpsi",accEffSubResultName);
3120
3121   Update();
3122
3123   return realSpectra;
3124 }
3125
3126 //_____________________________________________________________________________
3127 AliAnalysisMuMuSpectra* AliAnalysisMuMu::ComputeYield(const char* type, const char* flavour)
3128 {
3129   if (!SIM())
3130   {
3131     AliError("Cannot compute corrected yield without associated MC file !");
3132     return 0x0;
3133   }
3134
3135   const char* accEffSubResultName="COUNTJPSI:1";
3136
3137   AliAnalysisMuMuSpectra* realSpectra = GetSpectra(type,flavour);
3138
3139   if ( !realSpectra )
3140   {
3141     AliError("could not get real spectra");
3142     return 0x0;
3143   }
3144
3145   if (!realSpectra->HasValue("CoffNofJpsi"))
3146   {
3147     if (!CorrectSpectra(type,flavour))
3148     {
3149       AliError("Could not get corrected spectra");
3150       return 0x0;
3151     }
3152   }
3153
3154   AliAnalysisMuMuSpectra* simSpectra = SIM()->GetSpectra(type,flavour);
3155
3156   if ( !simSpectra)
3157   {
3158     AliErrorClass("could not get sim spectra");
3159     return 0x0;
3160   }
3161
3162   Double_t nofCMUL7 = CC()->GetSum(Form("trigger:CMUL7-B-NOPF-MUON/event:PSALL"));
3163   Double_t nofCINT7 = CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD/event:PSALL"));
3164   Double_t nofCINT7w0MUL = CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD&0MUL/event:PSALL"));
3165
3166   AliAnalysisMuMuBinning* binning = realSpectra->Binning();
3167   TObjArray* bins = binning->CreateBinObjArray();
3168   TIter nextBin(bins);
3169   AliAnalysisMuMuBinning::Range* bin;
3170   Int_t i(0);
3171   AliAnalysisMuMuJpsiResult* r;
3172
3173   while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
3174   {
3175     r = static_cast<AliAnalysisMuMuJpsiResult*>(realSpectra->BinContentArray()->At(i));
3176
3177     StdoutToAliDebug(1,std::cout << "bin=";r->Print(););
3178
3179     AliAnalysisMuMuJpsiResult* rsim = static_cast<AliAnalysisMuMuJpsiResult*>(simSpectra->BinContentArray()->At(i));
3180
3181     Double_t mbeq = nofCINT7w0MUL / ( nofCINT7 * nofCMUL7);
3182     Double_t mbeqError = mbeq * AliAnalysisMuMuResult::ErrorABC( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
3183                                                                 nofCINT7,TMath::Sqrt(nofCINT7),
3184                                                                 nofCMUL7,TMath::Sqrt(nofCMUL7));
3185
3186     r->Set("Fnorm",nofCINT7/nofCINT7w0MUL,(nofCINT7/nofCINT7w0MUL)*AliAnalysisMuMuResult::ErrorAB( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
3187                                                                                                 nofCINT7,TMath::Sqrt(nofCINT7)));
3188
3189     Double_t yield =  r->GetValue("CorrNofJpsi") * mbeq;
3190
3191     Double_t yieldError = yield * AliAnalysisMuMuResult::ErrorAB( r->GetValue("CorrNofJpsi"), r->GetErrorStat("CorrNofJpsi"),
3192                                                                  mbeq,mbeqError);
3193
3194     r->Set("YJpsi",yield,yieldError);
3195
3196     r->Set("NofInputJpsi",rsim->GetValue("NofInputJpsi",accEffSubResultName),rsim->GetErrorStat("NofInputJpsi",accEffSubResultName));
3197     r->Set("AccEffJpsi",rsim->GetValue("AccEffJpsi",accEffSubResultName),rsim->GetErrorStat("AccEffJpsi",accEffSubResultName));
3198
3199     ++i;
3200   }
3201
3202   delete bins;
3203
3204   Update();
3205
3206   return realSpectra;
3207 }
3208
3209 //_____________________________________________________________________________
3210 AliAnalysisMuMuSpectra* AliAnalysisMuMu::RABy(const char* realFile, const char* simFile, const char* type,
3211                                                const char* direction)
3212 {
3213   /// Compute the RAB...
3214   Double_t rapidityShift = 0.465;// 0.5*TMath::Log(208.0/82.0);
3215   const Double_t sqrts=5.023;
3216   const Double_t ymax=TMath::Log(sqrts*1000.0/3.096916);
3217   const Double_t tab = 0.093e-6; // nb^-1
3218   const Double_t tabError = 0.0035E-6; // nb^-1
3219   const char* accEffSubResultName="COUNTJPSI:1";
3220
3221   TF1 ydist("ydist","[0]*TMath::Exp(-(x*x)/(2.0*0.39*0.39))",0.,0.5);
3222   ydist.SetParameter(0,1.);
3223
3224   //Normalization to the values presented by Zaida and Rosana on January 11th 2013 https://indico.cern.ch/conferenceDisplay.py?confId=224985 slide 22
3225   // Normalization is done in the rapidity range 2.75<y<3.25 where Rosanas values is 230.8+212.1
3226   Double_t y1_norma= 2.75/ymax;
3227   Double_t y2_norma= 3.25/ymax;
3228   Double_t normalization = 0.25*(230.8+212.1)/ydist.Integral(y1_norma, y2_norma);
3229   ydist.SetParameter(0,normalization);
3230 //  AliInfoClass(Form("ymax=%e normalization=%f",ymax,ydist.Integral(y1_norma, y2_norma)));
3231
3232   AliAnalysisMuMu real(realFile);
3233   AliAnalysisMuMu sim(simFile);
3234
3235
3236   AliAnalysisMuMuSpectra* realSpectra = static_cast<AliAnalysisMuMuSpectra*>(real.MC()->GetObject(Form("/PSALL/CMUL7-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-%s",type)));
3237   AliAnalysisMuMuSpectra* simSpectra = static_cast<AliAnalysisMuMuSpectra*>(sim.MC()->GetObject(Form("/ALL/CMULLO-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-%s",type)));
3238
3239   if ( !realSpectra )
3240   {
3241     AliErrorClass("could not get real spectra");
3242     return 0x0;
3243   }
3244
3245   if ( !simSpectra)
3246   {
3247     AliErrorClass("could not get sim spectra");
3248     return 0x0;
3249   }
3250
3251   AliAnalysisMuMuSpectra* corrSpectra = static_cast<AliAnalysisMuMuSpectra*>(realSpectra->Clone());
3252   corrSpectra->Correct(*simSpectra,"Jpsi",accEffSubResultName);
3253
3254   Double_t nofCMUL7 = real.CC()->GetSum(Form("trigger:CMUL7-B-NOPF-MUON/event:PSALL"));
3255   Double_t nofCINT7 = real.CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD/event:PSALL"));
3256   Double_t nofCINT7w0MUL = real.CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD&0MUL/event:PSALL"));
3257
3258   AliAnalysisMuMuBinning* binning = realSpectra->Binning();
3259   TObjArray* bins = binning->CreateBinObjArray();
3260   TIter nextBin(bins);
3261   AliAnalysisMuMuBinning::Range* bin;
3262   Int_t i(0);
3263   AliAnalysisMuMuJpsiResult* r;
3264
3265   Int_t n = bins->GetLast();
3266
3267   TObjArray finalBins(n+1);
3268   finalBins.SetOwner(kTRUE);
3269
3270   TObjArray finalResults(n+1);
3271   finalResults.SetOwner(kFALSE);
3272
3273   while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
3274   {
3275     Double_t ylowlab = bin->Xmin();
3276     Double_t yhighlab = bin->Xmax();
3277
3278     Double_t ylowcms, yhighcms;
3279     Double_t ylownorm, yhighnorm;
3280
3281     if ( bin->IsIntegrated() )
3282     {
3283       ylowlab = -4;
3284       yhighlab = -2.5;
3285     }
3286
3287     if ( strcmp(direction,"pPb")==0 )
3288     {
3289       ylowcms = TMath::Abs(yhighlab) -  rapidityShift;
3290       yhighcms = TMath::Abs(ylowlab) - rapidityShift;
3291       ylownorm = ylowcms/ymax;
3292       yhighnorm = yhighcms/ymax;
3293     }
3294     else
3295     {
3296       ylowcms = ylowlab - rapidityShift;
3297       yhighcms = yhighlab - rapidityShift;
3298       ylownorm = -yhighcms/ymax;
3299       yhighnorm = -ylowcms/ymax;
3300     }
3301
3302
3303     Double_t brsigmapp = ydist.Integral(ylownorm,yhighnorm);
3304     Double_t brsigmappError = 0.0; // FIXME
3305
3306     AliInfoClass(Form("y range : LAB %f ; %f CMS %f ; %f -> ynorm : %f ; %f -> BR x sigmapp = %f",
3307                       ylowlab,yhighlab,ylowcms,yhighcms,ylownorm,yhighnorm,brsigmapp));
3308
3309     r = static_cast<AliAnalysisMuMuJpsiResult*>(corrSpectra->BinContentArray()->At(i)->Clone());
3310
3311     AliAnalysisMuMuJpsiResult* rsim = static_cast<AliAnalysisMuMuJpsiResult*>(simSpectra->BinContentArray()->At(i));
3312
3313     Double_t mbeq = nofCINT7w0MUL / ( nofCINT7 * nofCMUL7);
3314     Double_t mbeqError = mbeq * AliAnalysisMuMuResult::ErrorABC( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
3315                                          nofCINT7,TMath::Sqrt(nofCINT7),
3316                                          nofCMUL7,TMath::Sqrt(nofCMUL7));
3317
3318     r->Set("Fnorm",nofCINT7/nofCINT7w0MUL,(nofCINT7/nofCINT7w0MUL)*AliAnalysisMuMuResult::ErrorAB( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
3319                                                                       nofCINT7,TMath::Sqrt(nofCINT7)));
3320
3321     Double_t yield =  r->GetValue("CorrNofJpsi") * mbeq;
3322
3323     Double_t yieldError = yield * AliAnalysisMuMuResult::ErrorAB( r->GetValue("CorrNofJpsi"), r->GetErrorStat("CorrNofJpsi"),
3324                                           mbeq,mbeqError);
3325
3326     r->Set(Form("Y%sJpsi",direction),yield,yieldError);
3327
3328     Double_t raa = yield/(tab*brsigmapp);
3329     Double_t raaError = AliAnalysisMuMuResult::ErrorABC(yield,yieldError,
3330                                                         tab,tabError,
3331                                                         brsigmapp,brsigmappError);
3332     r->Set(Form("R%sJpsi",direction),raa,raaError);
3333
3334     r->Set("NofInputJpsi",rsim->GetValue("NofInputJpsi",accEffSubResultName),rsim->GetErrorStat("NofInputJpsi",accEffSubResultName));
3335     r->Set("AccEffJpsi",rsim->GetValue("AccEffJpsi",accEffSubResultName),rsim->GetErrorStat("AccEffJpsi",accEffSubResultName));
3336
3337     AliAnalysisMuMuBinning::Range* bincm = new AliAnalysisMuMuBinning::Range(bin->What(),bin->Quantity(),ylowcms,yhighcms);
3338
3339     r->SetBin(*bincm);
3340
3341     finalBins.Add(bincm);
3342     finalResults.Add(r);
3343
3344     ++i;
3345   }
3346
3347   delete bins;
3348
3349   AliAnalysisMuMuSpectra* spectra = new AliAnalysisMuMuSpectra(type,direction);
3350
3351   for ( i = 0; i <= n; ++i )
3352   {
3353     Int_t j(i);
3354     if ( strcmp(direction,"pPb")==0 )
3355     {
3356       j = n-i;
3357     }
3358
3359     r = static_cast<AliAnalysisMuMuJpsiResult*>(finalResults.At(j));
3360
3361     bin = static_cast<AliAnalysisMuMuBinning::Range*>(finalBins.At(j));
3362
3363     spectra->AdoptResult(*bin,r);
3364   }
3365
3366
3367   delete corrSpectra;
3368
3369   return spectra;
3370 }
3371
3372 //_____________________________________________________________________________
3373 void AliAnalysisMuMu::ComputeJpsi(const char* runlist, const char* prefix, const char* what, const char* binningFlavour)
3374 {
3375   /// Call the Jpsi method for each file
3376
3377   std::vector<int> runs;
3378   AliAnalysisTriggerScalers::ReadIntegers(runlist,runs);
3379
3380   for ( std::vector<int>::size_type i = 0; i < runs.size(); ++i )
3381   {
3382     Int_t runNumber = runs[i];
3383
3384     TString filename(Form("RUNBYRUN/%s_%09d.saf.root",prefix,runNumber));
3385
3386     if ( gSystem->AccessPathName(filename.Data())==kTRUE ) continue;
3387
3388     AliAnalysisMuMu m(filename.Data());
3389     m.Jpsi(what,binningFlavour);
3390   }
3391 }
3392
3393 //_____________________________________________________________________________
3394 AliAnalysisMuMuSpectra*
3395 AliAnalysisMuMu::CombineSpectra(const char* runlist,
3396                                 const char* realPrefix,
3397                                 const char* simPrefix,
3398                                 const char* quantity,
3399                                 const char* variable,
3400                                 const char* binningFlavour)
3401 {
3402   std::vector<int> runs;
3403   AliAnalysisTriggerScalers::ReadIntegers(runlist,runs);
3404
3405   std::vector<double> vw;
3406   std::vector<AliAnalysisMuMuSpectra*> vspectra;
3407
3408   for ( std::vector<int>::size_type i = 0; i < runs.size(); ++i )
3409   {
3410     Int_t runNumber = runs[i];
3411
3412     TString filename(Form("RUNBYRUN/%s_%09d.saf.root",realPrefix,runNumber));
3413
3414     AliAnalysisMuMu mreal(filename.Data());
3415
3416     std::cout << filename.Data() << std::flush << std::endl;
3417
3418     if ( !mreal.CC() || !mreal.MC() )
3419     {
3420       AliErrorClass(Form("Have to skip %s CC=%p MC=%p",filename.Data(),mreal.CC(),mreal.MC()));
3421       continue;
3422     }
3423
3424     TGraph* g = static_cast<TGraph*>(mreal.MC()->GetObject("/FNORM/GRAPHS/NMBeqBest2"));
3425
3426     if (!g)
3427     {
3428       mreal.ComputeFnorm();
3429       g  = static_cast<TGraph*>(mreal.MC()->GetObject("/FNORM/GRAPHS/NMBeqBest2"));
3430     }
3431
3432     if (!g) continue;
3433
3434     if ( TMath::Nint(g->GetX()[0]) != runNumber )
3435     {
3436       AliErrorClass("Wrong run number in NMBeqBest2 graph !");
3437       continue;
3438     }
3439
3440     filename.Form("RUNBYRUN/%s_%09d.saf.root",simPrefix,runNumber);
3441
3442     AliAnalysisMuMu msim(filename.Data());
3443
3444     TString spectraName(Form("/ALL/CMULLO-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-%s",variable));
3445     if (strlen(binningFlavour))
3446     {
3447       spectraName += "-";
3448       spectraName += binningFlavour;
3449     }
3450
3451     AliAnalysisMuMuSpectra* s = static_cast<AliAnalysisMuMuSpectra*>(msim.MC()->GetObject(spectraName.Data()));
3452
3453     if (!s)
3454     {
3455       msim.Jpsi(quantity,binningFlavour);
3456
3457       s = static_cast<AliAnalysisMuMuSpectra*>(msim.MC()->GetObject(spectraName.Data()));
3458     }
3459
3460     if (!s)
3461     {
3462       AliErrorClass(Form("Could not find spectra for sim %s",filename.Data()));
3463       continue;
3464     }
3465
3466     vw.push_back(g->GetY()[0]);
3467     vspectra.push_back(static_cast<AliAnalysisMuMuSpectra*>(s->Clone()));
3468   }
3469
3470   Double_t sumw(0.0);
3471
3472   for ( std::vector<int>::size_type i = 0; i < vw.size(); ++i )
3473   {
3474     sumw += vw[i];
3475     vspectra[i]->Print();
3476   }
3477
3478   // normalize the weights
3479   TList list;
3480   list.SetOwner(kTRUE);
3481
3482   for ( std::vector<int>::size_type i = 0; i < vw.size(); ++i )
3483   {
3484     vw[i] /= sumw;
3485
3486     vspectra[i]->SetWeight(vw[i]);
3487
3488     if ( i > 0 ) list.Add(vspectra[i]);
3489   }
3490
3491   std::cout << "before merging" << std::endl;
3492   for ( std::vector<int>::size_type i = 0; i < vw.size(); ++i )
3493   {
3494     std::cout << " ---------------- " << i << std::endl;
3495
3496     vspectra[i]->Print();
3497   }
3498
3499   vspectra[0]->Merge(&list);
3500
3501
3502   return vspectra[0];
3503 }
3504
3505
3506 //_____________________________________________________________________________
3507 TGraph* AliAnalysisMuMu::ResultEvolution(const char* runlist, const char* realPrefix, const char* simPrefix, const char* what, Bool_t forceRecomputation)
3508 {
3509   std::vector<int> runs;
3510   AliAnalysisTriggerScalers::ReadIntegers(runlist,runs);
3511
3512   TGraphErrors* g = new TGraphErrors(runs.size());
3513
3514   TString direction("Pbp");
3515   TString check(realPrefix);
3516   check.ToUpper();
3517
3518   if (check.Contains("LHC13B") ||
3519     check.Contains("LHC13C") ||
3520       check.Contains("LHC13D") ||
3521       check.Contains("LHC13E")
3522       )
3523   {
3524     direction = "pPb";
3525   }
3526
3527   AliInfoClass(Form("direction %s",direction.Data()));
3528
3529   Double_t mean(0.0);
3530   Double_t v1(0.0);
3531
3532   TString subResultName("");
3533
3534   TString swhat(what);
3535
3536   TObjArray* parts = swhat.Tokenize(":");
3537
3538   if ( parts->GetEntries() > 1 )
3539   {
3540     subResultName = static_cast<TObjString*>(parts->At(1))->String();
3541   }
3542
3543   delete parts;
3544
3545   for ( std::vector<int>::size_type i = 0; i < runs.size(); ++i )
3546   {
3547     Int_t runNumber = runs[i];
3548
3549     AliDebugClass(1,Form("RUN %09d",runNumber));
3550
3551     TString realFile(Form("RUNBYRUN/%s_%09d.saf.root",realPrefix,runNumber));
3552
3553     TString simFileName(Form("RUNBYRUN/%s_%09d.saf.root",simPrefix,runNumber));
3554
3555     TString simFile(simFileName);
3556
3557     TString resultName(Form("%s%sJpsi",what,direction.Data()));
3558
3559     if ( swhat == "Fnorm")
3560     {
3561       resultName = "Fnorm";
3562     }
3563     else if ( swhat.Contains("Acc") )
3564     {
3565       resultName.ReplaceAll(direction,"");
3566     }
3567
3568     AliAnalysisMuMu* mreal(0x0);
3569
3570     if ( !swhat.Contains("Acc"))
3571     {
3572         mreal = new AliAnalysisMuMu(realFile);
3573     }
3574
3575     AliAnalysisMuMuSpectra* real(0x0);
3576
3577     if ( mreal )
3578     {
3579       real = static_cast<AliAnalysisMuMuSpectra*>(mreal->MC()->GetObject("/PSALL/CMUL7-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-INTEGRATED"));
3580
3581       if (!real || forceRecomputation)
3582       {
3583         mreal->Jpsi();
3584         real = static_cast<AliAnalysisMuMuSpectra*>(mreal->MC()->GetObject("/PSALL/CMUL7-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-INTEGRATED"));
3585         if (!real)
3586         {
3587           AliErrorClass(Form("Could not get real spectra for run %d",runNumber));
3588           delete mreal;
3589           return 0x0;
3590         }
3591       }
3592     }
3593
3594     AliAnalysisMuMu msim(simFile);
3595
3596     AliAnalysisMuMuSpectra* sim = static_cast<AliAnalysisMuMuSpectra*>(msim.MC()->GetObject("/ALL/CMULLO-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-INTEGRATED"));
3597
3598     if (!sim || forceRecomputation)
3599     {
3600       msim.SetEventSelectionList("ALL");
3601       msim.Jpsi();
3602       sim = static_cast<AliAnalysisMuMuSpectra*>(msim.MC()->GetObject("/ALL/CMULLO-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-INTEGRATED"));
3603       if (!sim)
3604       {
3605         AliErrorClass(Form("Could not get sim spectra for run %d",runNumber));
3606         delete mreal;
3607         return 0x0;
3608       }
3609     }
3610
3611     AliAnalysisMuMuSpectra* corrected = 0x0;
3612
3613     if ( real )
3614     {
3615         corrected = AliAnalysisMuMu::RABy(realFile.Data(),simFile.Data(),"",direction.Data());
3616     }
3617     else
3618     {
3619       corrected = sim;
3620     }
3621
3622     AliAnalysisMuMuResult* result = static_cast<AliAnalysisMuMuResult*>(corrected->BinContentArray()->First());
3623
3624 //    result->Print();
3625
3626     Double_t value = result->GetValue(resultName.Data(),subResultName.Data());
3627     Double_t error = result->GetErrorStat(resultName.Data(),subResultName.Data());
3628
3629     g->SetPoint(i,runNumber,value);
3630     g->SetPointError(i,1,error);
3631
3632     Double_t n = 1.0;
3633
3634     if ( mreal )
3635     {
3636       n = mreal->CC()->GetSum(Form("trigger:CMUL7-B-NOPF-MUON/event:PSALL/run:%d",runNumber));
3637     }
3638
3639     Double_t w(0.0);
3640     if ( error > 0 ) w = 1.0/error/error;
3641
3642     mean += value*w;
3643     v1 += w;
3644
3645 //    std::cout << result->SubResults() << " " << result->GetError(resultName.Data()) << std::endl;
3646
3647     delete mreal;
3648   }
3649
3650   gStyle->SetOptFit(1);
3651   g->Draw("alp");
3652   g->Fit("pol0");
3653   g->SetTitle("");
3654   g->GetXaxis()->SetTitle("Run number");
3655   g->GetXaxis()->SetNoExponent();
3656   if ( TString(what) ==  "Y" )
3657   {
3658     g->GetYaxis()->SetTitle("J/#psi yield");
3659   }
3660   else if ( TString(what) == "R" )
3661   {
3662     g->GetYaxis()->SetTitle(Form("R_{%s}^{J/#psi}",direction.Data()));
3663   }
3664   else if ( TString(what).Contains("Acc") )
3665   {
3666     g->GetYaxis()->SetTitle("Acc#timesEff_{J/#psi}");
3667   }
3668   else if ( TString(what).Contains("Fnorm") )
3669   {
3670     g->GetYaxis()->SetTitle("CINT7 / CINT7&0MUL");
3671   }
3672
3673   if (CompactGraphs())
3674   {
3675     AliAnalysisMuMuGraphUtil::Compact(*g);
3676   }
3677
3678   mean /= v1;
3679
3680   AliInfoClass(Form("Weighted Mean %e",mean));
3681
3682   return g;
3683 }
3684