PWG/muondep/AliAnalysisMuMuFnorm.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 /// AliAnalysisMuMuFnorm : class to encapsulate computation(s)
  18 /// of the normalisation factor used to get the equivalent
  19 /// number of MB events from the number of CMUL triggers
  20 ///
  21 /// The computed objects are stored within a AliMergeableCollection
  22 /// with 3 subdirectories, dependinf on their type
  23 ///
  24 /// /GRAPHS/
  25 /// /RESULTS/
  26 /// /HISTOS/
  27 ///
  28 /// author: Laurent Aphecetche (Subatech)
  29
  30 #include "AliAnalysisMuMuFnorm.h"
  31
  32 #include "AliAnalysisMuMuGraphUtil.h"
  33 #include "AliAnalysisMuMuResult.h"
  34 #include "AliAnalysisTriggerScalers.h"
  35 #include "AliCounterCollection.h"
  36 #include "AliLog.h"
  37 #include "AliMergeableCollection.h"
  38 #include "Riostream.h"
  39 #include "TAxis.h"
  40 #include "TCanvas.h"
  41 #include "TGraphErrors.h"
  42 #include "TH1F.h"
  43 #include "TList.h"
  44 #include "TMap.h"
  45 #include "TMath.h"
  46 #include "TObjArray.h"
  47 #include "TObjString.h"
  48 #include "TPaveText.h"
  49 #include "TStyle.h"
  50 #include <cassert>
  51 #include <numeric>
  52
  53 ClassImp(AliAnalysisMuMuFnorm)
  54
  55 //_____________________________________________________________________________
  56 AliAnalysisMuMuFnorm::AliAnalysisMuMuFnorm(AliCounterCollection& cc,
  57                                            AliAnalysisMuMuFnorm::ETriggerType refTriggerType,
  58                                            const char* ocdbpath,
  59                                            Bool_t compactGraphs) :
  60 TObject(),
  61 fCounterCollection(cc),
  62 fMergeableCollection(0x0),
  63 fIsOwner(kTRUE),
  64 fOCDBPath(ocdbpath),
  65 fIsCompactGraphs(compactGraphs),
  66 fReferenceTriggerType(refTriggerType)
  67 {
  68   // ctor
  69
  70
  71 }
  72
  73 //_____________________________________________________________________________
  74 AliAnalysisMuMuFnorm::~AliAnalysisMuMuFnorm()
  75 {
  76   // dtor
  77   if ( fIsOwner )
  78   {
  79     delete fMergeableCollection;
  80   }
  81 }
  82
  83 //_____________________________________________________________________________
  84 void AliAnalysisMuMuFnorm::ComputeFnorm()
  85 {
  86   /// Compute the CMUL to CINT ratio(s)
  87   ///
  88   /// Using offline method
  89   ///   - in one go CINT/CMUL
  90   ///   - in two steps CINT/CMSL and CMSL/CMUL
  91   ///
  92   /// Using scaler method
  93   ///   - bare scaler values
  94   ///   - scaler values corrected for pile-up
  95   ///   - scaler values corrected for pile-up and physics selection
  96
  97   const ETriggerType triggerTypes[] = { kMB, kMUL, kMSL, kMSH };
  98   const Bool_t trueFalse[] = { kTRUE, kFALSE };
  99
 100   for ( Int_t i = 0; i < 4; ++i )
 101   {
 102     for ( Int_t pileup = 0; pileup < 2; ++pileup )
 103     {
 104       for ( Int_t ps = 0; ps < 2; ++ps )
 105       {
 106         ComputeNofEvents(triggerTypes[i],trueFalse[pileup],ps);
 107       }
 108     }
 109   }
 110
 111   ComputeFnormOffline(1,kFALSE,0);
 112   ComputeFnormOffline(1,kFALSE,1);
 113   ComputeFnormOffline(1,kTRUE,1);
 114
 115   ComputeFnormOffline(2,kFALSE,0);
 116   ComputeFnormOffline(2,kFALSE,1);
 117   ComputeFnormOffline(2,kTRUE,1);
 118
 119 //  ComputeFnormOffline(2,kFALSE,2);
 120 //  ComputeFnormOffline(2,kTRUE,2);
 121
 122   ComputeFnormScalers(kFALSE,0);
 123   ComputeFnormScalers(kTRUE,0);
 124   ComputeFnormScalers(kTRUE,1);
 125 //  ComputeFnormScalers(kTRUE,2);
 126
 127   WeightedMeanGraphs("Offline");
 128   WeightedMeanGraphs("Scalers");
 129   WeightedMeanGraphs("FnormOffline2PUPS,FnormOffline1PUPS","FnormOffline12PUPS");
 130
 131   WeightedMeanGraphs("FnormOffline2PUPS,FnormScalersPUPS","FnormBest2");
 132
 133   ComputeGraphRelDif("FnormOffline2PUPS","FnormScalersPUPS");
 134
 135   ComputeGraphRelDif("FnormOffline2PUPS","FnormOffline2");
 136   ComputeGraphRelDif("FnormOffline2PUPS","FnormOffline2PS");
 137
 138   ComputeGraphRelDif("CorrectionPSMB","CorrectionPSMUL");
 139
 140 //  for ( Int_t i = 0; i < 4; ++i )
 141 ///  {
 142     TString triggerEvents;
 143
 144 //  triggerEvents.Form("NofEvent%sPUPS",GetTriggerTypeName(triggerTypes[i]).Data());
 145   triggerEvents.Form("NofEvent%sPUPS",GetTriggerTypeName(fReferenceTriggerType).Data());
 146
 147   MultiplyGraphs(triggerEvents.Data(),"FnormBest2","NMBeqBest2");
 148
 149     MultiplyGraphs(triggerEvents.Data(),"FnormOffline2PUPS","NMBeqOffline2PUPS");
 150   MultiplyGraphs(triggerEvents.Data(),"FnormScalersPUPS","NMBeqScalersPUPS");
 151 //  }
 152
 153 //  MultiplyGraphs(Form("NofEvent%sPUTS",GetTriggerTypeName(fReferenceTriggerType).Data()),"FnormOffline2PUTS","NMBeqOffline2PUTS");
 154 //  MultiplyGraphs(Form("NofEvent%sPUTS",GetTriggerTypeName(fReferenceTriggerType).Data()),"FnormOffline2TS","NMBeqOffline2TS");
 155
 156   ComputeResultsFromGraphs();
 157
 158   AliAnalysisMuMuResult* result = GetResult("Fnorm");
 159   if (result)
 160   {
 161     result->Exclude("*");
 162     result->Include("FnormBest2");
 163   }
 164 }
 165
 166 //_____________________________________________________________________________
 167 void AliAnalysisMuMuFnorm::ComputeCorrectionFactors(Int_t eventSelectionCorrected)
 168 {
 169   /// Compute individual graphs for the correction factors (PS_CMUL, PS_CINT,
 170   /// F_pile-up,PS_CINT/PS_CMUL) used in the computation of (some) Fnorm factors
 171   ///
 172
 173   TString graphName(Form("CorrectionGlobal%s",GetEventSelectionName(eventSelectionCorrected).Data()));;
 174
 175   if ( GetGraph(graphName) )
 176   {
 177     // insure we compute it only once
 178     return;
 179   }
 180
 181   AliDebug(2,"");
 182
 183   std::vector<double> vx;
 184   std::vector<double> vxerr;
 185   std::vector<double> vy;
 186   std::vector<double> vyerr;
 187
 188   std::vector<double> vyGlobal;
 189   std::vector<double> vyGlobalErr;
 190
 191   const ETriggerType triggerTypes[] = { kMB, kMUL, kMSL, kMSH };
 192
 193   for ( Int_t i = 0; i < 4; ++i )
 194   {
 195     ComputeEventSelectionGraph(triggerTypes[i],eventSelectionCorrected);
 196     ComputePileUpGraph(triggerTypes[i],eventSelectionCorrected);
 197   }
 198
 199   TGraphErrors* gPSCINT = GetGraph(Form("Correction%s%s",GetEventSelectionName(eventSelectionCorrected).Data(),GetTriggerTypeName(AliAnalysisMuMuFnorm::kMB).Data()));
 200
 201   TGraphErrors* gPSCMUL = GetGraph(Form("Correction%s%s", GetEventSelectionName(eventSelectionCorrected).Data(), GetTriggerTypeName(AliAnalysisMuMuFnorm::kMUL).Data()));
 202
 203   TGraphErrors* gPU = GetGraph(Form("CorrectionPU%s%s",GetEventSelectionName(eventSelectionCorrected).Data(), GetTriggerTypeName(AliAnalysisMuMuFnorm::kMB).Data()));
 204
 205   if ( !gPSCINT || !gPSCMUL || !gPU )
 206   {
 207     AliError("Did not get the relevant graphs. Cannot work");
 208     return;
 209   }
 210
 211   for ( Int_t i = 0; i < gPSCINT->GetN(); ++i )
 212   {
 213     Double_t x,y,yerr,yGlobal,yGlobalErr;
 214
 215     gPSCINT->GetPoint(i,x,y);
 216
 217     if ( fIsCompactGraphs )
 218     {
 219       x = TString(gPSCINT->GetXaxis()->GetBinLabel(i)).Atoi();
 220     }
 221
 222     yGlobal = gPSCINT->GetY()[i] * gPU->GetY()[i] / gPSCMUL->GetY()[i];
 223
 224     yGlobalErr = yGlobal*AliAnalysisMuMuResult::ErrorABC(gPSCINT->GetY()[i],gPSCINT->GetEY()[i],
 225                                                          gPSCMUL->GetY()[i],gPSCMUL->GetEY()[i],
 226                                                          gPU->GetY()[i],gPU->GetEY()[i]);
 227
 228     y = gPSCINT->GetY()[i] / gPSCMUL->GetY()[i];
 229     yerr = y * AliAnalysisMuMuResult::ErrorAB(gPSCINT->GetY()[i],gPSCINT->GetEY()[i],
 230                                               gPSCMUL->GetY()[i],gPSCMUL->GetEY()[i]);
 231
 232     vx.push_back(x);
 233     vxerr.push_back(gPSCINT->GetEX()[i]);
 234
 235     vyGlobal.push_back(yGlobal);
 236     vyGlobalErr.push_back(yGlobalErr);
 237
 238     vy.push_back(y);
 239     vyerr.push_back(yerr);
 240   }
 241
 242   TString name(Form("Correction%sRatio",GetEventSelectionName(eventSelectionCorrected).Data()));
 243   TString title(Form("%s_MB/%s_%s",GetEventSelectionName(eventSelectionCorrected).Data(),GetTriggerTypeName(fReferenceTriggerType).Data(),
 244                      GetEventSelectionName(eventSelectionCorrected).Data()));
 245
 246   CreateAndAddGraph(name,title,vx,vxerr,vy,vyerr);
 247
 248   title = TString::Format("%s_MB x Fpile-up / %s_%s ",GetEventSelectionName(eventSelectionCorrected).Data(),GetTriggerTypeName(fReferenceTriggerType).Data(),GetEventSelectionName(eventSelectionCorrected).Data());
 249
 250   CreateAndAddGraph(graphName,title,vx,vxerr,vyGlobal,vyGlobalErr);
 251 }
 252
 253 //_____________________________________________________________________________
 254 void AliAnalysisMuMuFnorm::ComputeFnormOffline(Int_t nstep, Bool_t pileUpCorrected, Int_t eventSelectionCorrected)
 255 {
 256   /// Compute MB to CMUL ratio using offline method, either in 1 or 2 steps
 257
 258   TString name("FnormOffline");
 259   TString title("Computed using offline information");
 260
 261   if (nstep==1)
 262   {
 263     name += "1";
 264     title += " in one step (CINT/CINT&0MUL)";
 265   }
 266   else
 267   {
 268     name += "2";
 269     title += " in two steps (CMSL/CMSL&0MUL x CINT/CINT&MSL)";
 270   }
 271
 272   if (pileUpCorrected)
 273   {
 274     name += "PU";
 275     title += " with pile-up correction";
 276   }
 277   if (eventSelectionCorrected==1)
 278   {
 279     name += "PS";
 280     title += " with (ps) purity corrections";
 281   }
 282   else if ( eventSelectionCorrected==2 )
 283   {
 284     name += "TS";
 285     title += " with (ts) purity corrections";
 286   }
 287
 288   if ( GetGraph(name) )
 289   {
 290     // insure we're computing it only once
 291     return;
 292   }
 293
 294   AliDebug(2,name);
 295
 296   std::vector<double> vx;
 297   std::vector<double> vxerr;
 298   std::vector<double> vy;
 299   std::vector<double> vyerr;
 300
 301   const std::set<int>& runs = RunNumbers();
 302
 303   for ( std::set<int>::const_iterator it = runs.begin(); it != runs.end(); ++it )
 304   {
 305     Int_t runNumber = *it;
 306
 307     TString mbTrigger = GetTriggerClassName(kMB,runNumber);
 308     TString muonTrigger = GetTriggerClassName(kMSL,runNumber);
 309     TString dimuonTrigger = GetTriggerClassName(kMUL,runNumber);
 310
 311     if (!mbTrigger.Length())
 312     {
 313       AliError(Form("Cannot get MB trigger for run %d",runNumber));
 314       continue;
 315     }
 316
 317     Double_t nofMB = GetSum(mbTrigger.Data(),runNumber,eventSelectionCorrected);
 318     Double_t nofMSL(0.0);
 319     Double_t nofMSLw0MUL(0.0);
 320
 321     if ( nstep==2 )
 322     {
 323       nofMSL = GetSum(muonTrigger.Data(),runNumber,eventSelectionCorrected);
 324       nofMSLw0MUL = GetSum(Form("%s&0MUL",muonTrigger.Data()),runNumber,eventSelectionCorrected);
 325     }
 326
 327     Double_t nofMBw0MUL = GetSum(Form("%s&0MUL",mbTrigger.Data()),runNumber,eventSelectionCorrected);
 328     Double_t nofMBw0MSL = GetSum(Form("%s&0MSL",mbTrigger.Data()),runNumber,eventSelectionCorrected);
 329
 330     if ( !nofMBw0MUL ) continue;
 331     if ( !nofMBw0MSL && nstep == 2 ) continue;
 332
 333     Double_t purityMB(1.0);
 334     Double_t purityMBerror(0.0);
 335
 336     if ( eventSelectionCorrected > 0 )
 337     {
 338       ComputeEventSelectionGraph(kMB,eventSelectionCorrected);
 339
 340       TGraphErrors* gps = GetGraph(Form("Correction%s%s",GetEventSelectionName(eventSelectionCorrected).Data(),GetTriggerTypeName(kMB).Data()));
 341
 342       GetValueAndErrorFromGraph(gps,runNumber,purityMB,purityMBerror);
 343     }
 344
 345     Double_t pileUpFactor(1.0);
 346     Double_t pileUpFactorError(0.0);
 347
 348     if (pileUpCorrected)
 349     {
 350       ComputePileUpGraph(kMB,eventSelectionCorrected);
 351
 352       TGraphErrors* gpu = GetGraph(Form("CorrectionPU%s%s",GetEventSelectionName(eventSelectionCorrected).Data(),GetTriggerTypeName(kMB).Data()));
 353
 354       GetValueAndErrorFromGraph(gpu,runNumber,pileUpFactor,pileUpFactorError);
 355
 356       nofMB *= pileUpFactor;
 357     }
 358
 359     double value = nofMBw0MUL > 0.0 ? nofMB/nofMBw0MUL : 0.0;
 360     double error = value*AliAnalysisMuMuResult::ErrorABC(nofMB,TMath::Sqrt(nofMB),
 361                                                               nofMBw0MUL,TMath::Sqrt(nofMBw0MUL),
 362                                                               pileUpFactor,pileUpFactorError);
 363
 364     if ( nstep == 2 )
 365     {
 366       value = (nofMB/nofMSLw0MUL)*(nofMSL/nofMBw0MSL);
 367
 368       if ( runNumber == 196310 )
 369       {
 370         AliDebug(1,Form("RUN %09d %d-%d-%d value=%e nofMB %e nofMSLw0MUL %e nofMSL %e nofMBw0MSL %e",
 371                         runNumber,nstep,pileUpCorrected,eventSelectionCorrected,
 372                         value,nofMB,nofMSLw0MUL,nofMSL,nofMBw0MSL));
 373       }
 374
 375       error = value*AliAnalysisMuMuResult::ErrorABCD(nofMB,TMath::Sqrt(nofMB),
 376                                                           nofMSLw0MUL,TMath::Sqrt(nofMSLw0MUL),
 377                                                           nofMSL,TMath::Sqrt(nofMSL),
 378                                                           nofMBw0MSL,TMath::Sqrt(nofMBw0MSL));
 379     }
 380
 381     if ( value > 0.0 )
 382     {
 383       vx.push_back(1.0*runNumber);
 384       vxerr.push_back(0.5);
 385       vy.push_back(value);
 386       vyerr.push_back(error);
 387     }
 388   }
 389
 390
 391   CreateAndAddGraph(name,title,vx,vxerr,vy,vyerr);
 392 }
 393
 394
 395 //_____________________________________________________________________________
 396 void AliAnalysisMuMuFnorm::ComputeFnormScalers(Bool_t pileUpCorrected,
 397                                                Int_t eventSelectionCorrected)
 398 {
 399   /// Compute the MB to CMUL ratio using the scalers method (from OCDB)
 400   ///
 401   /// i.e. Fnorm = L0B(MB) x PS(MB) x Fpile-up / ( L0B(REF) x PS(REF) )
 402   ///
 403   /// where MB is the minbias trigger
 404   /// REF is the fReferenceTrigger
 405   /// and PS is the fraction of events selected by the physics selection
 406   ///
 407   /// The correction factor (the two PS and one Fpile-up) are
 408   /// taken from graphs computed in other methods
 409   ///
 410
 411   TString name("FnormScalers");
 412   TString title("Computed using OCDB scalers");
 413
 414   if (pileUpCorrected)
 415   {
 416     name += "PU";
 417     title += " with pile-up correction";
 418   }
 419   if (eventSelectionCorrected==1)
 420   {
 421     name += "PS";
 422     title += " with (ps) purity corrections";
 423   }
 424   if (eventSelectionCorrected==2)
 425   {
 426     name += "TS";
 427     title += " with (ts) purity corrections";
 428   }
 429
 430   if ( GetGraph(name) )
 431   {
 432     // insure we're computing it only once
 433     return;
 434   }
 435
 436   AliDebug(2,name);
 437
 438   // insure we have all the graphs we need to work
 439   ComputeTriggerL0B(kMB);
 440   ComputeTriggerL0B(fReferenceTriggerType);
 441
 442   const std::set<int>& runs = RunNumbers();
 443
 444   std::vector<double> vx;
 445   std::vector<double> vxerr;
 446   std::vector<double> vy;
 447   std::vector<double> vyerr;
 448
 449   Double_t purityREF(1.0);
 450   Double_t purityMB(1.0);
 451   Double_t purityREFerror(00);
 452   Double_t purityMBerror(0.0);
 453
 454   // compute the per run values
 455   for ( std::set<int>::const_iterator it = runs.begin(); it != runs.end(); ++it )
 456   {
 457     Int_t runNumber = *it;
 458
 459     TString mbTrigger = GetTriggerClassName(kMB,runNumber);
 460     TString refTrigger = GetTriggerClassName(fReferenceTriggerType,runNumber);
 461
 462     purityMB=purityREF=1.0;
 463     purityMBerror=purityREFerror=0.0;
 464
 465     if (eventSelectionCorrected>0)
 466     {
 467       ComputeEventSelectionGraph(kMB,eventSelectionCorrected);
 468       ComputeEventSelectionGraph(fReferenceTriggerType,eventSelectionCorrected);
 469
 470       TGraphErrors* gpsMB  = GetGraph(Form("Correction%s%s",GetEventSelectionName(eventSelectionCorrected).Data(),GetTriggerTypeName(kMB).Data()));
 471       TGraphErrors* gpsREF = GetGraph(Form("Correction%s%s",GetEventSelectionName(eventSelectionCorrected).Data(),GetTriggerTypeName(fReferenceTriggerType).Data()));
 472
 473       GetValueAndErrorFromGraph(gpsMB,runNumber,purityMB,purityMBerror);
 474       GetValueAndErrorFromGraph(gpsREF,runNumber,purityREF,purityREFerror);
 475     }
 476
 477     if (purityMB<=0.0)
 478     {
 479       AliError(Form("Got purity=%e for MB for run %9d",purityMB,runNumber));
 480       continue;
 481     }
 482
 483     TGraphErrors* gl0bMB = GetGraph(Form("L0B%s",GetTriggerTypeName(kMB).Data()));
 484     TGraphErrors* gl0bREF = GetGraph(Form("L0B%s",GetTriggerTypeName(fReferenceTriggerType).Data()));
 485
 486     Double_t L0bMB,L0bMBError;
 487     Double_t L0bREF,L0bREFError;
 488
 489     GetValueAndErrorFromGraph(gl0bMB,runNumber,L0bMB,L0bMBError);
 490     GetValueAndErrorFromGraph(gl0bREF,runNumber,L0bREF,L0bREFError);
 491
 492     Double_t pileUpFactor(1.0);
 493     Double_t pileUpFactorError(0.0);
 494
 495     if (pileUpCorrected)
 496     {
 497       TGraphErrors* gpu = GetGraph((Form("CorrectionPU%s%s",GetEventSelectionName(eventSelectionCorrected).Data(),GetTriggerTypeName(kMB).Data())));
 498
 499       GetValueAndErrorFromGraph(gpu,runNumber,pileUpFactor,pileUpFactorError);
 500     }
 501
 502     Double_t value;
 503     Double_t error;
 504
 505     ScalerFnorm(value,error,
 506                 L0bREF,purityREF,purityREFerror,
 507                 L0bMB,purityMB,purityMBerror,
 508                 pileUpFactor,pileUpFactorError);
 509
 510     if ( value > 0.0 )
 511     {
 512       vx.push_back(1.0*runNumber);
 513       vxerr.push_back(0.5);
 514       vy.push_back(value);
 515       vyerr.push_back(error);
 516     }
 517   }
 518
 519   CreateAndAddGraph(name,title,vx,vxerr,vy,vyerr);
 520 }
 521
 522 //_____________________________________________________________________________
 523 void AliAnalysisMuMuFnorm::ComputeGraphRelDif(const char* a, const char* b) const
 524 {
 525   // compute dispersion of b versus a
 526   //
 527   // computed differences graphs are put into the GRAPHS/ directory
 528   // computed differences results are put into the HISTOS/ directory
 529
 530   TString name(Form("RelDif%svs%s",b,a));
 531
 532   if ( GetGraph(name) )
 533   {
 534     // insure we compute it only once
 535     return;
 536   }
 537
 538   AliDebug(2,name);
 539
 540   TGraphErrors* ga = static_cast<TGraphErrors*>(MC()->GetObject(Form("/GRAPHS/%s",a)));
 541   TGraphErrors* gb = static_cast<TGraphErrors*>(MC()->GetObject(Form("/GRAPHS/%s",b)));
 542
 543   if (!ga)
 544   {
 545     AliError(Form("Cannot get graph for %s",a));
 546     return;
 547   }
 548
 549   if (!gb)
 550   {
 551     AliError(Form("Cannot get graph for %s",b));
 552     return;
 553   }
 554
 555   if ( ga->GetN() != gb->GetN() )
 556   {
 557     AliError(Form("Cannot work with different number of points in the graphs : %d vs %d",
 558                   ga->GetN(),gb->GetN()));
 559     return;
 560   }
 561
 562   TString title(Form("%s-%s (RelDif,%%)",b,a));
 563
 564   std::vector<double> vx;
 565   std::vector<double> vxerr;
 566   std::vector<double> vy;
 567   std::vector<double> vyerr;
 568
 569   for ( Int_t i = 0; i < ga->GetN(); ++i )
 570   {
 571     Double_t xa,xb,ya,yb;
 572
 573     ga->GetPoint(i,xa,ya);
 574     gb->GetPoint(i,xb,yb);
 575
 576     if ( xa != xb )
 577     {
 578       AliError(Form("Incompatible graphs : got xa=%e and xb=%e",xa,xb));
 579       return;
 580     }
 581
 582     Double_t newvalue = 0.0;
 583
 584     if ( TMath::Abs(xa) > 1E-12 )
 585     {
 586       newvalue = 100.0*( yb - ya ) / ya;
 587     }
 588
 589     Double_t yerr = newvalue*AliAnalysisMuMuResult::ErrorAB(ya,ga->GetEY()[i],
 590                                                             yb,gb->GetEY()[i]);
 591
 592     if ( fIsCompactGraphs )
 593     {
 594       xa = TString(ga->GetXaxis()->GetBinLabel(i+1)).Atoi()*1.0;
 595     }
 596
 597     vx.push_back(xa);
 598     vxerr.push_back(0.5);
 599     vy.push_back(newvalue);
 600     vyerr.push_back(yerr);
 601
 602   }
 603
 604   CreateAndAddGraph(name,title,vx,vxerr,vy,vyerr);
 605
 606   // FIXME : fill here an histogram from the graph to get the
 607   // weight of 1/e2 ?
 608   //  h->Fill(newvalue,1.0/(yerr*yerr));
 609   //  MC()->Adopt("/HISTOS/",h);
 610
 611   //  AliAnalysisMuMuResult* r = GetRunIntegratedResult(*g,"FnormDispersion");
 612   //  if (r)
 613   //  {
 614   //    if (!dispersion)
 615   //    {
 616   //      dispersion = new AliAnalysisMuMuResult("FnormDispersion");
 617   //    }
 618   //    dispersion->AdoptSubResult(r);
 619   //    if ( !TString(g->GetName()).BeginsWith("Fnorm") )
 620   //    {
 621   //      dispersion->Exclude(r->Alias());
 622   //    }
 623   //  }
 624 }
 625
 626 //_____________________________________________________________________________
 627 void AliAnalysisMuMuFnorm::ComputePileUpGraph(ETriggerType tt, Int_t eventSelectionCorrected)
 628 {
 629   /// Compute the per-run graph of pile-up factor
 630
 631   TString graphName(Form("CorrectionPU%s%s",GetEventSelectionName(eventSelectionCorrected).Data(),GetTriggerTypeName(tt).Data()));
 632
 633   if ( GetGraph(graphName) )
 634   {
 635     // insure we're computing it only once
 636     return;
 637   }
 638
 639   AliDebug(2,graphName);
 640
 641   std::vector<double> vx;
 642   std::vector<double> vxerr;
 643   std::vector<double> vy;
 644   std::vector<double> vyerr;
 645
 646   const std::set<int>& runs = RunNumbers();
 647
 648   AliAnalysisTriggerScalers ts(runs,OCDBPath().Data());
 649
 650   for ( std::set<int>::const_iterator it = runs.begin(); it != runs.end(); ++it )
 651   {
 652     Int_t runNumber = *it;
 653
 654     Double_t pileUpFactor(1.0);
 655     Double_t pileUpFactorError(0.0);
 656     Double_t purity(1.0);
 657     Double_t purityError(0.0);
 658
 659     TString triggerClassName = GetTriggerClassName(tt,runNumber);
 660
 661     if ( triggerClassName.Length()==0 )
 662     {
 663       AliError(Form("Unknown trigger type %d",tt));
 664       return;
 665     }
 666
 667     if (eventSelectionCorrected)
 668     {
 669       GetPurity(triggerClassName.Data(),runNumber,purity,purityError,eventSelectionCorrected);
 670     }
 671     ts.GetPileUpFactor(runNumber,triggerClassName.Data(),purity,pileUpFactor,pileUpFactorError);
 672
 673     vx.push_back(runNumber);
 674     vxerr.push_back(0.5);
 675     vy.push_back(pileUpFactor);
 676     vyerr.push_back(pileUpFactorError);
 677   }
 678
 679   TString title(Form("Pile-up correction factor for trigger %s",GetTriggerTypeName(tt).Data()));
 680
 681   if (eventSelectionCorrected)
 682   {
 683     title += "( L0BRate corrected by event selection";
 684     title += GetEventSelectionName(eventSelectionCorrected);
 685     title += " accept fraction)";
 686   }
 687
 688   CreateAndAddGraph(graphName,title,vx,vxerr,vy,vyerr);
 689 }
 690
 691 //_____________________________________________________________________________
 692 void AliAnalysisMuMuFnorm::ComputeEventSelectionGraph(ETriggerType tt, Int_t eventSelectionCorrected)
 693 {
 694   /// Compute the per-run graph of physics selection accept fraction
 695   /// for the given trigger
 696
 697   TString graphName(Form("Correction%s%s",GetEventSelectionName(eventSelectionCorrected).Data(), GetTriggerTypeName(tt).Data()));
 698
 699   if (GetGraph(graphName))
 700   {
 701     // insure we're computing it only once
 702     return;
 703   }
 704
 705   AliDebug(2,graphName);
 706
 707   std::vector<double> vx;
 708   std::vector<double> vxerr;
 709   std::vector<double> vy;
 710   std::vector<double> vyerr;
 711
 712   const std::set<int>& runs = RunNumbers();
 713
 714   AliAnalysisTriggerScalers ts(runs,OCDBPath().Data());
 715
 716   for ( std::set<int>::const_iterator it = runs.begin(); it != runs.end(); ++it )
 717   {
 718     Int_t runNumber = *it;
 719
 720     Double_t purity, purityError;
 721
 722     TString triggerClassName = GetTriggerClassName(tt,runNumber);
 723
 724     if ( triggerClassName.Length()==0 )
 725     {
 726       AliError(Form("Unknown trigger type %d",tt));
 727       return;
 728     }
 729
 730     GetPurity(triggerClassName.Data(),runNumber,purity,purityError,eventSelectionCorrected);
 731
 732     vx.push_back(runNumber);
 733     vxerr.push_back(0.5);
 734     vy.push_back(purity);
 735     vyerr.push_back(purityError);
 736   }
 737
 738   TString title(Form("Fraction of events accepted by the event selection %s for trigger %s",GetTriggerTypeName(tt).Data(),
 739                      GetEventSelectionName(eventSelectionCorrected).Data()));
 740
 741   CreateAndAddGraph(graphName,title,vx,vxerr,vy,vyerr);
 742 }
 743
 744
 745 //_____________________________________________________________________________
 746 void AliAnalysisMuMuFnorm::ComputeResultsFromGraphs()
 747 {
 748   // Compute one single value for each graph, by weighting by the fraction
 749   // of events in each run
 750   // do this for certain groups of graphs
 751
 752   TObjArray groups;
 753   groups.SetOwner(kTRUE);
 754
 755   groups.Add(new TObjString("Fnorm"));
 756   groups.Add(new TObjString("NMBeq"));
 757   groups.Add(new TObjString("Correction"));
 758   groups.Add(new TObjString("RelDif"));
 759
 760   TList* objectList = MC()->CreateListOfObjectNames("/GRAPHS/");
 761
 762   TIter nextGroup(&groups);
 763   TObjString* grp;
 764
 765   TIter next(objectList);
 766   TObjString* str(0x0);
 767
 768   while ( ( grp = static_cast<TObjString*>(nextGroup()) ) )
 769   {
 770     TString oname(Form("/RESULTS/%s",grp->String().Data()));
 771
 772     if ( MC()->GetObject(oname) )
 773     {
 774       // delete if we have it already so we can replace it
 775       MC()->Remove(oname);
 776     }
 777
 778     AliAnalysisMuMuResult* result = new AliAnalysisMuMuResult(grp->String());
 779
 780     MC()->Adopt("/RESULTS/",result);
 781
 782     next.Reset();
 783
 784     while ( ( str = static_cast<TObjString*>(next()) ) )
 785     {
 786       if ( ! ( str->String().BeginsWith(grp->String() ) ) ) continue;
 787
 788       TGraphErrors* g = GetGraph(str->String());
 789
 790       if (!g) continue;
 791
 792       AliAnalysisMuMuResult* sub = GetRunIntegratedResult(*g);
 793
 794       if ( !sub )
 795       {
 796         AliError(Form("Could not get result for graph %s",g->GetName()));
 797       }
 798       if ( sub )
 799       {
 800         result->AdoptSubResult(sub);
 801       }
 802     }
 803   }
 804
 805   delete objectList;
 806 }
 807
 808 //_____________________________________________________________________________
 809 void AliAnalysisMuMuFnorm::ComputeNofEvents(ETriggerType triggerType,
 810                                             Bool_t pileUpCorrected,
 811                                             Int_t eventSelectionCorrected)
 812 {
 813   /// Compute trigger fractions
 814
 815   TString graphName(Form("NofEvent%s%s%s",GetTriggerTypeName(triggerType).Data(),
 816                          pileUpCorrected ? "PU" : "",
 817                          GetEventSelectionName(eventSelectionCorrected).Data()));
 818
 819   if ( GetGraph(graphName) )
 820   {
 821     // compute it only once
 822     return;
 823   }
 824
 825   ComputeCorrectionFactors(eventSelectionCorrected);
 826
 827   TString gpsname(Form("Correction%s%s",GetEventSelectionName(eventSelectionCorrected).Data(),GetTriggerTypeName(triggerType).Data()));
 828   TString gpuname(Form("CorrectionPU%s%s",GetEventSelectionName(eventSelectionCorrected).Data(), GetTriggerTypeName(triggerType).Data()));
 829
 830   TGraphErrors* gPS = GetGraph(gpsname);
 831
 832   if (!gPS)
 833   {
 834     AliError(Form("Could not find %s",gpsname.Data()));
 835     return;
 836   }
 837
 838   TGraphErrors* gPU = GetGraph(gpuname);
 839
 840   if (!gPU)
 841   {
 842     AliError(Form("Could not find %s",gpuname.Data()));
 843     return;
 844   }
 845
 846   AliDebug(2,graphName);
 847
 848   std::vector<double> vx;
 849   std::vector<double> vxerr;
 850   std::vector<double> vy;
 851   std::vector<double> vyerr;
 852
 853   const std::set<int>& runs = RunNumbers();
 854
 855   Int_t i(0);
 856
 857   for ( std::set<int>::const_iterator it = runs.begin(); it != runs.end(); ++it )
 858   {
 859     Int_t runNumber = *it;
 860
 861     TString triggerClassName = GetTriggerClassName(triggerType,runNumber);
 862
 863     if ( triggerClassName.Length() )
 864     {
 865       Double_t n = GetSum(triggerClassName,runNumber,0);
 866
 867       vx.push_back(runNumber);
 868       vxerr.push_back(0.5);
 869
 870       assert(runNumber==TMath::Nint(gPU->GetX()[i]));
 871
 872       Double_t y(n);
 873       Double_t y1(1.0);
 874       Double_t y2(1.0);
 875       Double_t e1(0);
 876       Double_t e2(0);
 877
 878       if ( pileUpCorrected )
 879       {
 880         y1 = gPU->GetY()[i];
 881         e1 = gPU->GetEY()[i];
 882       }
 883
 884       if ( eventSelectionCorrected > 0 )
 885       {
 886         y2 = gPS->GetY()[i];
 887         e2 = gPS->GetEY()[i];
 888       }
 889
 890       y *= y1*y2;
 891
 892       AliDebug(2,Form("RUN %09d n %e y1 %e y2 %e y% e",runNumber,n,y1,y2,y));
 893
 894       Double_t yerr = y*AliAnalysisMuMuResult::ErrorABC( n, TMath::Sqrt(n),
 895                                                         y1, e1,
 896                                                         y2, e2);
 897
 898       vy.push_back(y);
 899       vyerr.push_back(yerr);
 900
 901       ++i;
 902     }
 903   }
 904
 905   TString title(Form("Number of event of trigger %s",GetTriggerTypeName(triggerType).Data()));
 906
 907   CreateAndAddGraph(graphName,title,vx,vxerr,vy,vyerr);
 908 }
 909
 910 //_____________________________________________________________________________
 911 void AliAnalysisMuMuFnorm::ComputeTriggerFractions(ETriggerType triggerType,
 912                                                    Bool_t physicsSelectionCorrected)
 913 {
 914   /// Compute trigger fractions
 915
 916   TString graphName(Form("Fractions%s%s",GetTriggerTypeName(triggerType).Data(),physicsSelectionCorrected ? "PS" : ""));
 917
 918   if ( GetGraph(graphName) )
 919   {
 920     // compute it only once
 921     return;
 922   }
 923
 924   AliDebug(2,graphName);
 925
 926   std::vector<double> vx;
 927   std::vector<double> vxerr;
 928   std::vector<double> vy;
 929   std::vector<double> vyerr;
 930
 931   const std::set<int>& runs = RunNumbers();
 932
 933   Double_t n(0.0);
 934
 935   for ( std::set<int>::const_iterator it = runs.begin(); it != runs.end(); ++it )
 936   {
 937     Int_t runNumber = *it;
 938
 939     TString triggerClassName = GetTriggerClassName(triggerType,runNumber);
 940
 941     if ( triggerClassName.Length() )
 942     {
 943       n += GetSum(triggerClassName,runNumber,physicsSelectionCorrected);
 944     }
 945   }
 946
 947   if ( n <= 0.0 )
 948   {
 949     AliWarning(Form("Got zero trigger for %s",GetTriggerTypeName(triggerType).Data()));
 950     return;
 951   }
 952
 953   for ( std::set<int>::const_iterator it = runs.begin(); it != runs.end(); ++it )
 954   {
 955     Int_t runNumber = *it;
 956
 957     TString triggerClassName = GetTriggerClassName(triggerType,runNumber);
 958
 959     vx.push_back(runNumber);
 960     vxerr.push_back(0.5);
 961     Double_t y = GetSum(triggerClassName,runNumber,physicsSelectionCorrected);
 962     vy.push_back(y/n);
 963     vyerr.push_back( (y/n) * AliAnalysisMuMuResult::ErrorAB( y,TMath::Sqrt(y),
 964                                                         n, TMath::Sqrt(n)));
 965   }
 966
 967
 968   TString title(Form("Fraction of event of trigger %s",GetTriggerTypeName(triggerType).Data()));
 969
 970   CreateAndAddGraph(graphName,title,vx,vxerr,vy,vyerr);
 971 }
 972
 973 //_____________________________________________________________________________
 974 void AliAnalysisMuMuFnorm::ComputeTriggerL0B(ETriggerType triggerType)
 975 {
 976   /// Compute trigger L0B
 977
 978   std::vector<double> vx;
 979   std::vector<double> vxerr;
 980   std::vector<double> vy;
 981   std::vector<double> vyerr;
 982
 983   TString graphName(Form("L0B%s",GetTriggerTypeName(triggerType).Data()));
 984
 985   if ( GetGraph(graphName) )
 986   {
 987     // insure we're computing it only once
 988     return;
 989   }
 990
 991   AliDebug(2,graphName);
 992
 993   const std::set<int>& runs = RunNumbers();
 994
 995   AliAnalysisTriggerScalers ts(runs,OCDBPath().Data());
 996
 997   for ( std::set<int>::const_iterator it = runs.begin(); it != runs.end(); ++it )
 998   {
 999     Int_t runNumber = *it;
1000
1001     TString triggerClassName = GetTriggerClassName(triggerType,runNumber);
1002
1003     AliAnalysisTriggerScalerItem* item = ts.GetTriggerScaler(runNumber,"L0B",triggerClassName);
1004     if (!item) continue;
1005
1006     vx.push_back(runNumber);
1007     vxerr.push_back(0.5);
1008
1009     Double_t y = item->Value();
1010
1011     vy.push_back(y);
1012
1013     vyerr.push_back( TMath::Sqrt(y) );
1014   }
1015
1016   TString title(Form("L0B of trigger %s",GetTriggerTypeName(triggerType).Data()));
1017
1018   CreateAndAddGraph(graphName,title,vx,vxerr,vy,vyerr);
1019 }
1020
1021 ////_____________________________________________________________________________
1022 //void AliAnalysisMuMuFnorm::ComputeTSGraph(ETriggerType tt)
1023 //{
1024 //  /// Compute the per-run graph of physics selection accept fraction x track
1025 //  /// accept fraction for the given trigger
1026 //
1027 //  TString graphName(Form("CorrectionTS%s",GetTriggerTypeName(tt).Data()));
1028 //
1029 //  if (GetGraph(graphName))
1030 //  {
1031 //    // insure we're computing it only once
1032 //    return;
1033 //  }
1034 //
1035 //  AliDebug(1,graphName);
1036 //
1037 //  std::vector<double> vx;
1038 //  std::vector<double> vxerr;
1039 //  std::vector<double> vy;
1040 //  std::vector<double> vyerr;
1041 //
1042 //  const std::set<int>& runs = RunNumbers();
1043 //
1044 //  AliAnalysisTriggerScalers ts(runs,OCDBPath().Data());
1045 //
1046 //  for ( std::set<int>::const_iterator it = runs.begin(); it != runs.end(); ++it )
1047 //  {
1048 //    Int_t runNumber = *it;
1049 //
1050 //    Double_t purity, purityError;
1051 //
1052 //    TString triggerClassName = GetTriggerClassName(tt,runNumber);
1053 //
1054 //    if ( triggerClassName.Length()==0 )
1055 //    {
1056 //      AliError(Form("Unknown trigger type %d",tt));
1057 //      return;
1058 //    }
1059 //
1060 //    GetPurity(triggerClassName.Data(),runNumber,purity,purityError,"OFFLINE1");
1061 //
1062 //    vx.push_back(runNumber);
1063 //    vxerr.push_back(0.5);
1064 //    vy.push_back(purity);
1065 //    vyerr.push_back(purityError);
1066 //  }
1067 //
1068 //  TString title(Form("Fraction of events accepted by the physics selection x track selection for trigger %s",GetTriggerTypeName(tt).Data()));
1069 //
1070 //  CreateAndAddGraph(graphName,title,vx,vxerr,vy,vyerr);
1071 //}
1072 //
1073
1074 //_____________________________________________________________________________
1075 TGraphErrors* AliAnalysisMuMuFnorm::CreateAndAddGraph(const TString& name,
1076                                                       const TString& title,
1077                                                       const std::vector<double>& vx,
1078                                                       const std::vector<double>& vxerr,
1079                                                       const std::vector<double>& vy,
1080                                                       const std::vector<double>& vyerr) const
1081 {
1082   /// Creates a graph and adds it to our mergeable collection
1083
1084   TGraphErrors* g = new TGraphErrors(vx.size(),&vx[0],&vy[0],&vxerr[0],&vyerr[0]);
1085   g->SetName(name.Data());
1086   g->SetTitle(title.Data());
1087
1088   if  (fIsCompactGraphs)
1089   {
1090     AliAnalysisMuMuGraphUtil::Compact(*g);
1091   }
1092
1093   g->GetXaxis()->SetNoExponent();
1094 //  g->GetXaxis()->SetTitle("Run number");
1095
1096   TPaveText* text = new TPaveText(0.70,0.70,0.89,0.89,"NDC");
1097   text->SetBorderSize(0);
1098   text->SetFillColor(0);
1099   text->AddText(Form("Mean %e",g->GetMean(2)));
1100   text->AddText(Form("RMS  %e",g->GetRMS(2)));
1101   g->GetListOfFunctions()->Add(text);
1102
1103   MC()->Adopt("/GRAPHS/",g);
1104   return g;
1105 }
1106
1107 //_____________________________________________________________________________
1108 AliMergeableCollection* AliAnalysisMuMuFnorm::DetachMC()
1109 {
1110   // let go the ownership of our mergeable collection
1111   fIsOwner = kFALSE;
1112   return fMergeableCollection;
1113 }
1114
1115 //_____________________________________________________________________________
1116 void AliAnalysisMuMuFnorm::DrawWith2Scales(const char* graphName1, const char* graphName2)
1117 {
1118   TGraphErrors* g1 = static_cast<TGraphErrors*>(GetGraph(graphName1)->Clone());
1119   TGraphErrors* g2 = static_cast<TGraphErrors*>(GetGraph(graphName2)->Clone());
1120
1121   if ( g1 && g2 )
1122   {
1123     gStyle->SetOptTitle(0);
1124
1125     AliAnalysisMuMuGraphUtil gu;
1126
1127     TCanvas* c = gu.DrawWith2Scales(*g1,*g2);
1128     c->Draw();
1129
1130     for ( Int_t i = 0; i < 2; ++i )
1131     {
1132       c->cd(i);
1133       gPad->SetLeftMargin(0.15);
1134       gPad->SetRightMargin(0.15);
1135     }
1136
1137     c->Update();
1138
1139     //    TGraphErrors* g = new TGraphErrors(g1->GetN());
1140     //
1141     //    Double_t check(0.0);
1142     //
1143     //    for ( Int_t i = 0; i < g->GetN(); ++i )
1144     //    {
1145     //      Double_t y = g1->GetY()[i]*g2->GetY()[i];
1146     //
1147     //      check += g2->GetY()[i];
1148     //
1149     //      g->SetPoint(i,g2->GetX()[i],y);
1150     //      g->SetPointError(i,g2->GetEX()[i],
1151     //                       y*AliAnalysisMuMuResult::ErrorAB(g1->GetY()[i],g1->GetEY()[i],
1152     //                                                        g2->GetY()[i],g2->GetEY()[i]));
1153     //    }
1154     //
1155     //    new TCanvas;
1156     //
1157     //    g->Draw("ap");
1158     //
1159     //    AliInfo(Form("check: %e g mean %e rms %e",check,g->GetMean(2),g->GetRMS(2)));
1160
1161     /*
1162
1163     // g1 vs g2
1164
1165     TGraphErrors* g = new TGraphErrors(g1->GetN());
1166
1167     for ( Int_t i = 0; i < g->GetN(); ++i )
1168     {
1169       g->SetPoint(i,g2->GetY()[i],g1->GetY()[i]);
1170       g->SetPointError(i,g2->GetEY()[i],g1->GetEY()[i]);
1171     }
1172
1173     new TCanvas;
1174
1175     g->Draw("ap");
1176
1177     AliInfo(Form("g mean %e rms %e",g->GetMean(2),g->GetRMS(2)));
1178
1179      */
1180
1181   }
1182
1183 }
1184
1185 //_____________________________________________________________________________
1186 TString AliAnalysisMuMuFnorm::GetEventSelectionName(Int_t eventSelectionCorrected) const
1187 {
1188   if ( eventSelectionCorrected == 1 )
1189   {
1190     return "PS";
1191   }
1192   else if ( eventSelectionCorrected == 2 )
1193   {
1194     return "TS";
1195   }
1196   return "";
1197 }
1198
1199 //_____________________________________________________________________________
1200 TGraphErrors* AliAnalysisMuMuFnorm::GetGraph(const char* name) const
1201 {
1202   // shortcut method to give access to one graph
1203
1204   TObject* o = MC()->GetObject(Form("/GRAPHS/%s",name));
1205
1206   if (!o) return 0x0;
1207
1208   if ( o->IsA() != TGraphErrors::Class() )
1209   {
1210     AliError(Form("Object %s is not of the expected type",o->GetName()));
1211     return 0x0;
1212   }
1213
1214   return static_cast<TGraphErrors*>(o);
1215 }
1216
1217 //_____________________________________________________________________________
1218 void AliAnalysisMuMuFnorm::GetPurity(const char* triggerClassName, Int_t runNumber,
1219                                      Double_t& value, Double_t& error, Int_t eventSelectionCorrected) const
1220 {
1221   /// Get the physics selection accept fraction for a given trigger
1222
1223   value=error=0.0;
1224
1225   TString runCondition;
1226
1227   if (runNumber>0)
1228   {
1229     runCondition.Form("/run:%d",runNumber);
1230   }
1231
1232   Double_t nall = fCounterCollection.GetSum(Form("trigger:%s/event:ALL%s",triggerClassName,runCondition.Data()));
1233
1234   TString ename;
1235
1236   if ( eventSelectionCorrected == 1 )
1237   {
1238     ename = "PSALL";
1239   }
1240   else if ( eventSelectionCorrected == 2 )
1241   {
1242     ename = "OFFLINE1";
1243   }
1244   else
1245   {
1246     value = 1.0;
1247     return;
1248
1249   }
1250   Double_t nps = fCounterCollection.GetSum(Form("trigger:%s/event:%s%s",
1251                                                 triggerClassName,ename.Data(),runCondition.Data()));
1252
1253   if ( nall <= 0.0 ) return;
1254
1255   value = nps/nall;
1256
1257   error = AliAnalysisMuMuResult::ErrorAB(nall,TMath::Sqrt(nall),nps,TMath::Sqrt(nps));
1258 }
1259
1260 //_____________________________________________________________________________
1261 AliAnalysisMuMuResult* AliAnalysisMuMuFnorm::GetResult(const char* name) const
1262 {
1263   // shortcut method to give access to one result
1264
1265   TObject* o = MC()->GetObject(Form("/RESULTS/%s",name));
1266
1267   if (!o) return 0x0;
1268
1269   if ( o->IsA() != AliAnalysisMuMuResult::Class() )
1270   {
1271     AliError(Form("Object %s is not of the expected type",o->GetName()));
1272     return 0x0;
1273   }
1274
1275   return static_cast<AliAnalysisMuMuResult*>(o);
1276 }
1277
1278 //______________________________________________________________________________
1279 AliAnalysisMuMuResult* AliAnalysisMuMuFnorm::GetRunIntegratedResult(const TGraphErrors& g, const char* basename)
1280 {
1281   /// get one value +- error for this graph (weighting the run-by-run values
1282   /// by the fraction of the number of triggers in that run)
1283   /// also the rms is computed.
1284
1285   Bool_t physicsSelectionCorrected(kFALSE);
1286
1287   if ( TString(g.GetName()).Contains("PS") )
1288   {
1289     physicsSelectionCorrected=kTRUE;
1290   }
1291
1292   ComputeTriggerFractions(fReferenceTriggerType,physicsSelectionCorrected);
1293
1294   TString fname(Form("Fractions%s%s",GetTriggerTypeName(fReferenceTriggerType).Data(),physicsSelectionCorrected ? "PS" : ""));
1295
1296   TGraphErrors* gTriggerFractions = GetGraph(fname);
1297
1298   StdoutToAliDebug(2,std::cout << g.GetName() << std::endl; g.Print(););
1299
1300   if (!gTriggerFractions)
1301   {
1302     AliError(Form("Did not find graph for %s",fname.Data()));
1303     return 0x0;
1304   }
1305
1306   Double_t mean = g.GetY()[0];
1307   Int_t n = g.GetN();
1308
1309   Double_t errorOnMean = g.GetEY()[0];
1310   Double_t rms = 0.0;
1311
1312   if ( n > 1 )
1313   {
1314     mean = errorOnMean = 0.0;
1315     Double_t d(0.0);
1316     Double_t v1(0.0);
1317     Double_t v2(0.0);
1318
1319     for ( Int_t i = 0; i < n; ++i )
1320     {
1321       Double_t y = g.GetY()[i];
1322
1323       Double_t weight = gTriggerFractions->GetY()[i]; // sum of weight should be 1.0
1324
1325       AliDebug(2,Form("%s i %3d y %e weight %e",g.GetName(),i,y,weight));
1326
1327       mean += y * weight;
1328
1329       v1 += weight;
1330       v2 += weight*weight;
1331     }
1332
1333     mean /= v1;
1334
1335     for ( Int_t i = 0; i < n; ++i )
1336     {
1337       Double_t weight = gTriggerFractions->GetY()[i]; // sum of weight should be 1.0
1338       Double_t y = g.GetY()[i];
1339
1340       d += (y-mean)*(y-mean)*weight;
1341     }
1342
1343     AliDebug(2,Form("v1=%e v2=%e d=%e",v1,v2,d));
1344
1345     if ( v1 <= 0) return 0x0;
1346
1347     errorOnMean = TMath::Sqrt((1.0/v1)*(1.0/(n-1))*d);
1348
1349     rms = TMath::Sqrt( (v1/(v1*v1-v2))*d );
1350   }
1351
1352   AliAnalysisMuMuResult* result = new AliAnalysisMuMuResult(g.GetName(),g.GetTitle());
1353
1354   result->Set(basename,mean,errorOnMean,rms);
1355
1356   if ( TString(g.GetName()) == "FnormScalersPUPS" )
1357   {
1358     AliDebug(1,Form("mean %e errorOnMean %e rms %e",mean,errorOnMean,rms));
1359     StdoutToAliDebug(1,result->Print("full"));
1360   }
1361
1362   return result;
1363 }
1364
1365 //_____________________________________________________________________________
1366 Double_t AliAnalysisMuMuFnorm::GetSum(const char* triggerClassName,
1367                                       Int_t runNumber,
1368                                       Int_t eventSelectionCorrected) const
1369 {
1370   TString condition(Form("trigger:%s/run:%d",triggerClassName,runNumber));
1371
1372   if (eventSelectionCorrected==1)
1373   {
1374     condition += "/event:PSALL";
1375   }
1376   else if ( eventSelectionCorrected == 2 )
1377   {
1378     condition += "/event:OFFLINE1";
1379   }
1380   else
1381   {
1382     condition += "/event:ALL";
1383   }
1384
1385   Double_t n = fCounterCollection.GetSum(condition.Data());
1386
1387   if (n<=0)
1388   {
1389     AliError(Form("Got no count for %s for run %d (physicsSelected:%d)",triggerClassName,runNumber,eventSelectionCorrected));
1390     return 0;
1391   }
1392
1393   return n;
1394 }
1395
1396 //_____________________________________________________________________________
1397 TString AliAnalysisMuMuFnorm::GetTriggerClassName(ETriggerType tt, Int_t runNumber) const
1398 {
1399   // get the triggerclass to for a given trigger type and run number
1400
1401   if ( tt == kMB )
1402   {
1403     return MBTriggerClassName(runNumber);
1404   }
1405   else if ( tt == kMUL )
1406   {
1407     return MULTriggerClassName(runNumber);
1408   }
1409   else if ( tt == kMSL)
1410   {
1411     return MSLTriggerClassName(runNumber);
1412   }
1413   else if ( tt == kMSH)
1414   {
1415     return MSHTriggerClassName(runNumber);
1416   }
1417   return "";
1418 }
1419
1420 //_____________________________________________________________________________
1421 TString AliAnalysisMuMuFnorm::GetTriggerTypeName(ETriggerType tt) const
1422 {
1423   // get the name of the trigger type
1424   if ( tt == kMB )
1425   {
1426     return "MB";
1427   }
1428   else if ( tt == kMUL )
1429   {
1430     return "MUL";
1431   }
1432   else if ( tt == kMSL)
1433   {
1434     return "MSL";
1435   }
1436   else if ( tt == kMSH)
1437   {
1438     return "MSH";
1439   }
1440   return "";
1441 }
1442
1443 //_____________________________________________________________________________
1444 void AliAnalysisMuMuFnorm::GetValueAndErrorFromGraph(TGraphErrors* graph,
1445                                                      Int_t runNumber,
1446                                                      Double_t& value,
1447                                                      Double_t& error) const
1448 {
1449   /// get (value,error) corresponding to run=runNumber.
1450   /// Works both for compact and non-compact graphs
1451
1452   value = TMath::Limits<Double_t>::Max();
1453   error = 0.0;
1454
1455   if (!graph) return;
1456
1457   TAxis* axis = graph->GetXaxis();
1458
1459   for ( Int_t i = 0; i < graph->GetN(); ++i )
1460   {
1461     Int_t rbin = TMath::Nint(graph->GetX()[i]);
1462     Int_t rlabel = TString(axis->GetBinLabel(i+1)).Atoi();
1463     if ( rbin == runNumber || rlabel == runNumber )
1464     {
1465       value = graph->GetY()[i];
1466       error = graph->GetEY()[i];
1467     }
1468   }
1469 }
1470
1471 //_____________________________________________________________________________
1472 AliMergeableCollection* AliAnalysisMuMuFnorm::MC() const
1473 {
1474   // get our mergeable collection
1475   if (!fMergeableCollection)
1476   {
1477     fMergeableCollection = new AliMergeableCollection("Fnorm",Form("MB to %s trigger normalization results",GetTriggerTypeName(fReferenceTriggerType).Data()));
1478   }
1479   return fMergeableCollection;
1480 }
1481
1482 //_____________________________________________________________________________
1483 TString AliAnalysisMuMuFnorm::MBTriggerClassName(Int_t runNumber) const
1484 {
1485   /// FIXME : find a better way ?
1486
1487   if ( TriggerClassnameTest("CINT7-B-NOPF-ALLNOTRD",runNumber) )
1488   {
1489     return "CINT7-B-NOPF-ALLNOTRD";
1490   }
1491   return "";
1492 }
1493
1494 //_____________________________________________________________________________
1495 TString AliAnalysisMuMuFnorm::MSHTriggerClassName(Int_t runNumber) const
1496 {
1497   /// FIXME : find a better way ?
1498
1499   if ( TriggerClassnameTest("CMSH7-B-NOPF-ALLNOTRD",runNumber) )
1500   {
1501     return "CMSH7-B-NOPF-ALLNOTRD";
1502   }
1503   else if ( TriggerClassnameTest("CMSH7-B-NOPF-MUON",runNumber) )
1504   {
1505     return "CMSH7-B-NOPF-MUON";
1506   }
1507   return "";
1508 }
1509
1510 //_____________________________________________________________________________
1511 TString AliAnalysisMuMuFnorm::MSLTriggerClassName(Int_t runNumber) const
1512 {
1513   /// FIXME : find a better way ?
1514
1515   if ( TriggerClassnameTest("CMSL7-B-NOPF-MUON",runNumber) )
1516   {
1517       return "CMSL7-B-NOPF-MUON";
1518   }
1519 //  else
1520 //    if ( TriggerClassnameTest("CMSL7-B-NOPF-ALLNOTRD",runNumber) )
1521 //  {
1522 //    return "CMSL7-B-NOPF-ALLNOTRD";
1523 //  }
1524   return "";
1525 }
1526
1527 //_____________________________________________________________________________
1528 void AliAnalysisMuMuFnorm::MultiplyGraphs(const char* g1name, const char* g2name, const char* name)
1529 {
1530   /// Make a new graph = g1*g2
1531   std::vector<double> vx;
1532   std::vector<double> vy;
1533   std::vector<double> vxerr;
1534   std::vector<double> vyerr;
1535
1536   TGraphErrors* g1 = GetGraph(g1name);
1537   TGraphErrors* g2 = GetGraph(g2name);
1538
1539   if (!g1)
1540   {
1541     AliError(Form("Could not get graph %s",g1name));
1542     return;
1543   }
1544
1545   if (!g2)
1546   {
1547     AliError(Form("Could not get graph %s",g2name));
1548     return;
1549   }
1550
1551   if ( g1->GetN() != g2->GetN() )
1552   {
1553     AliError(Form("Could not multiply incompatible graphs %d pts vs %d pts",g1->GetN(),g2->GetN()));
1554     return;
1555   }
1556
1557   for ( Int_t i = 0; i < g1->GetN(); ++i )
1558   {
1559     if ( g1->GetX()[i] != g2->GetX()[i] )
1560     {
1561       AliError(Form("Incompatible bin %d : %e vs %e",i,g1->GetX()[i],g2->GetX()[i]));
1562       return;
1563     }
1564
1565     vx.push_back(g1->GetX()[i]);
1566     vxerr.push_back(g1->GetEX()[i]);
1567
1568     Double_t y = g1->GetY()[i]*g2->GetY()[i];
1569     Double_t yerr = y*AliAnalysisMuMuResult::ErrorAB( g1->GetY()[i], g1->GetEY()[i],
1570                                                      g2->GetY()[i], g2->GetEY()[i]);
1571
1572     vy.push_back(y);
1573     vyerr.push_back(yerr);
1574   }
1575
1576   TString gname(name);
1577
1578   if ( gname.Length() == 0 )
1579   {
1580     gname = g1->GetName();
1581     gname += "x";
1582     gname += g2->GetName();
1583   }
1584
1585   TString title(Form("Product of %s by %s",g1->GetName(),g2->GetName()));
1586
1587   CreateAndAddGraph(gname,title,vx,vxerr,vy,vyerr);
1588 }
1589
1590 //_____________________________________________________________________________
1591 TString AliAnalysisMuMuFnorm::MULTriggerClassName(Int_t runNumber) const
1592 {
1593   /// FIXME : find a better way ?
1594
1595   if ( TriggerClassnameTest("CMUL7-B-NOPF-ALLNOTRD",runNumber) )
1596   {
1597     return "CMUL7-B-NOPF-ALLNOTRD";
1598   }
1599   else if ( TriggerClassnameTest("CMUL7-B-NOPF-MUON",runNumber) )
1600   {
1601     return "CMUL7-B-NOPF-MUON";
1602   }
1603   return "";
1604
1605 }
1606
1607 //_____________________________________________________________________________
1608 void AliAnalysisMuMuFnorm::Print(Option_t* opt) const
1609 {
1610   if ( fMergeableCollection )
1611   {
1612     fMergeableCollection->Print(opt);
1613   }
1614 }
1615
1616 //_____________________________________________________________________________
1617 std::set<int> AliAnalysisMuMuFnorm::RunNumbers() const
1618 {
1619   // Extract the run numbers from our counter collection
1620
1621   std::set<int> runset;
1622
1623   TString sruns = fCounterCollection.GetKeyWords("run");
1624   TObjArray* runs = sruns.Tokenize(",");
1625
1626   TIter next(runs);
1627   TObjString* s;
1628
1629   while ( ( s = static_cast<TObjString*>(next())) )
1630   {
1631     runset.insert(s->String().Atoi());
1632   }
1633   delete runs;
1634
1635   return runset;
1636 }
1637
1638 //_____________________________________________________________________________
1639 void AliAnalysisMuMuFnorm::ScalerFnorm(Double_t& value, Double_t& error,
1640                                        Double_t L0bREF, Double_t purityREF, Double_t purityREFerror,
1641                                        Double_t L0bMB, Double_t purityMB, double_t purityMBerror,
1642                                        Double_t pileUpFactor, Double_t pileUpFactorError)
1643 {
1644   /// Compute the MB to CMUL ratio and its associated error
1645
1646   value = error = 0.0;
1647
1648   value = L0bREF*purityREF;
1649
1650   if (!value) return;
1651
1652   value = L0bMB*purityMB*pileUpFactor/value;
1653
1654   error = value*AliAnalysisMuMuResult::ErrorABCDE(L0bREF,TMath::Sqrt(L0bREF),
1655                                                   purityREF,purityREFerror,
1656                                                   L0bMB,TMath::Sqrt(L0bMB),
1657                                                   purityMB,purityMBerror,
1658                                                   pileUpFactor,pileUpFactorError);
1659 }
1660
1661 //_____________________________________________________________________________
1662 void AliAnalysisMuMuFnorm::ShowFnorm(const TObjArray& a) const
1663 {
1664   /// Print and plot Fnorm values
1665   TIter next(&a);
1666   TGraphErrors* g;
1667
1668   while ( ( g = static_cast<TGraphErrors*>(next()) ) )
1669   {
1670     g->SetTitle(Form("Fnorm %s",g->GetTitle()));
1671
1672     new TCanvas(Form("c%s",g->GetName()));
1673
1674     if (fIsCompactGraphs)
1675     {
1676       AliAnalysisMuMuGraphUtil::Compact(*g);
1677     }
1678     else
1679     {
1680       g->GetXaxis()->SetNoExponent();
1681     }
1682     g->Draw("lpa");
1683
1684     Double_t y(0.0);
1685     Double_t yerr(0.0);
1686
1687     for ( Int_t i = 0; i < g->GetN(); ++i )
1688     {
1689       y += g->GetY()[i];
1690       Double_t e = ( y != 0.0  ? g->GetEY()[i]/y : 0.0);
1691
1692       yerr += e*e;
1693     }
1694
1695     y /= (g->GetN());
1696     yerr = TMath::Sqrt(yerr)*y;
1697
1698     AliInfo(Form("%30s graph %e +- %e (%5.2f %%) RMS %e",g->GetName(),
1699                  y,yerr,yerr*100/y,g->GetRMS(2)));
1700
1701   }
1702 }
1703
1704 //_____________________________________________________________________________
1705 Bool_t AliAnalysisMuMuFnorm::TriggerClassnameTest(const char* triggerClassName, Int_t runNumber) const
1706 {
1707   /// Check if we have counts for that trigger,run combination
1708
1709   TString runs = fCounterCollection.GetKeyWords("run");
1710
1711   if ( !runs.Contains(Form("%d",runNumber)) ) return kFALSE;
1712
1713   TString triggers = fCounterCollection.GetKeyWords("trigger");
1714
1715   if (!triggers.Contains(triggerClassName)) return kFALSE;
1716
1717   Double_t n = fCounterCollection.GetSum(Form("trigger:%s/run:%d",triggerClassName,runNumber));
1718
1719   return ( n > 0.0 );
1720 }
1721
1722 //_____________________________________________________________________________
1723 void
1724 AliAnalysisMuMuFnorm::WeightedMeanGraphs(const char* patternOrList, const char* graphName)
1725 {
1726   /// Sum the graphs which name matches pattern
1727   /// Sum is made using a weighted mean (each element is weighted by the inverse
1728   /// of its error squared)
1729
1730   TString spattern(patternOrList);
1731   TObjArray* slist(0x0);
1732
1733   if ( spattern.CountChar(',') )
1734   {
1735     // it's not a pattern but a list...
1736     slist = spattern.Tokenize(",");
1737     spattern = "";
1738   }
1739
1740   TList* objectList = MC()->CreateListOfObjectNames("/GRAPHS/");
1741   TIter next(objectList);
1742   TObjString* str(0x0);
1743   TObjArray selected;
1744   selected.SetOwner(kFALSE);
1745
1746   while ( ( str = static_cast<TObjString*>(next()) ) )
1747   {
1748     TGraphErrors* g = GetGraph(str->String());
1749
1750     if (!g) continue;
1751
1752     TString name(g->GetName());
1753
1754     if ( spattern.Length() >0 && !name.Contains(spattern.Data()) ) continue;
1755
1756     if ( slist && !slist->FindObject(name)) continue;
1757
1758     AliDebug(2,Form("Selected for sum : %s",name.Data()));
1759
1760     selected.Add(g);
1761   }
1762
1763   delete slist;
1764   delete objectList;
1765
1766   if (selected.GetLast()<0) return;
1767
1768   std::vector<double> vx;
1769   std::vector<double> vy;
1770   std::vector<double> vxerr;
1771   std::vector<double> vyerr;
1772
1773   Int_t npts = static_cast<TGraphErrors*>(selected.First())->GetN();
1774
1775   for ( Int_t ipoint = 0; ipoint < npts; ++ipoint )
1776   {
1777     Double_t x,xref,xerr;
1778     Double_t sum(0.0);
1779     Double_t sume2(0.0);
1780
1781     for ( Int_t igraph = 0; igraph <= selected.GetLast(); ++igraph )
1782     {
1783       TGraphErrors* g = static_cast<TGraphErrors*>(selected.At(igraph));
1784
1785       if ( g->GetN() != npts )
1786       {
1787         AliError(Form("Graph %s does not have the expected %d points",g->GetName(),npts));
1788         continue;
1789       }
1790       Double_t runNumber;
1791
1792       if ( fIsCompactGraphs )
1793       {
1794         runNumber = TString(g->GetXaxis()->GetBinLabel(ipoint+1)).Atoi()*1.0;
1795       }
1796       else
1797       {
1798         runNumber = g->GetX()[ipoint];
1799       }
1800
1801       if ( igraph == 0 )
1802       {
1803         xref = g->GetX()[ipoint];
1804         x = runNumber;
1805         xerr = g->GetEX()[ipoint];
1806       }
1807       else
1808       {
1809         if ( xref != g->GetX()[ipoint] )
1810         {
1811           AliError(Form("Cannot sum graphs with different axis : get %e and expected %e : %s vs %s",xref,x,selected.At(0)->GetName(),g->GetName()));
1812           return;
1813         }
1814       }
1815
1816       Double_t e2 = g->GetEY()[ipoint]*g->GetEY()[ipoint];
1817
1818       if ( e2 != 0.0 )
1819       {
1820         sum += g->GetY()[ipoint]/e2;
1821         sume2 += 1.0/e2;
1822       }
1823     }
1824
1825     if (sume2 != 0.0)
1826     {
1827       vx.push_back(x);
1828       vxerr.push_back(xerr);
1829       vy.push_back(sum/sume2);
1830       vyerr.push_back(TMath::Sqrt(1/sume2));
1831     }
1832   }
1833
1834   Int_t n = selected.GetEntries();
1835
1836   TString name(graphName);
1837   TString title(Form("Weighted mean from %d individual graphs",n));
1838
1839   if ( strlen(graphName) == 0)
1840   {
1841     name = TString::Format("WeightMeanFnorm%s",patternOrList);
1842     name.ReplaceAll(",","_");
1843     title = TString::Format("WeightMeanFnorm%s from %d individual graphs",patternOrList,n);
1844   }
1845
1846   CreateAndAddGraph(name,title,vx,vxerr,vy,vyerr);
1847 }
1848
1849