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