PWG/muondep/AliAnalysisTriggerScalers.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 // $Id$
  17
  18 //
  19 // AliAnalysisTriggerScalers : utility class to play with OCDB scalers.
  20 //
  21 // Can use it to e.g. :
  22 //
  23 // - get the integrated luminosity for a given period
  24 //   (see IntegratedLuminosity method)
  25 //
  26 // - plot the trigger rates (see PlotTriggerEvolution)
  27 //
  28 // Please note that this class is doing an OCDB "scan" (for as many run
  29 // as you put in the SetRunList method), so it can be really long
  30 // if you're using the raw:// OCDB. If you're planning on using this
  31 // class for anything more than a brief test, you'd better think
  32 // or making a local copy of the OCDB objects required by this class :
  33 //
  34 // GRP/CTP/Config
  35 // GRP/CTP/Scalers
  36 //
  37 // and also (to get the fill and period ranges, mainly for drawing)
  38 //
  39 // GRP/GRP/Data
  40 // GRP/GRP/LHCData
  41 //
  42 //
  43 // Please note also that this class evolved from a set a quick macros to
  44 // follow the luminosity and trigger rates during data taking to this stage.
  45 //
  46 // Now (Feb 2013) would probably be a good time to regroup a bit, think about it and
  47 // make it more general, more robust and just less "had-hoc"...
  48 //
  49 // Ideas for improvement :
  50 //
  51 // - make sure the "what" parameter mean the same thing in all methods
  52 //   and so can take the same values in all methods
  53 //
  54 // - get the lumi trigger name and cross-section from e.g. OADB instead of
  55 //   hard-coding them ?
  56 //
  57 // - organize the class so that the CTP Scalers and Config can be fetched
  58 //   from elsewhere (e.g. from a run based container available in some
  59 //   distant future in the ESDs/AODs ?)
  60 //
  61 // - robustify the PAC fraction computation
  62 //
  63 //
  64 // L. Aphecetche (Subatech)
  65 //
  66
  67 #include "AliAnalysisTriggerScalers.h"
  68 #include "AliCDBEntry.h"
  69 #include "AliCDBManager.h"
  70 #include "AliGRPObject.h"
  71 #include "AliLog.h"
  72 #include "AliTimeStamp.h"
  73 #include "AliTriggerBCMask.h"
  74 #include "AliTriggerClass.h"
  75 #include "AliTriggerCluster.h"
  76 #include "AliTriggerConfiguration.h"
  77 #include "AliTriggerDescriptor.h"
  78 #include "AliTriggerInput.h"
  79 #include "AliTriggerRunScalers.h"
  80 #include "AliTriggerScalers.h"
  81 #include "AliTriggerScalersESD.h"
  82 #include "AliTriggerScalersRecord.h"
  83 #include "Riostream.h"
  84 #include "TDatime.h"
  85 #include "TError.h"
  86 #include "TGraph.h"
  87 #include "TGraphErrors.h"
  88 #include "TH2.h"
  89 #include "TObjArray.h"
  90 #include "TObjString.h"
  91 #include <set>
  92 #include <map>
  93 #include "TStyle.h"
  94 #include "AliLHCData.h"
  95 #include "TMath.h"
  96 #include "TAxis.h"
  97 #include "TBox.h"
  98 #include "TCanvas.h"
  99 #include "TText.h"
 100 using std::make_pair;
 101 ClassImp(AliAnalysisTriggerScalers)
 102
 103 namespace {
 104
 105   //______________________________________________________________________________
 106   UChar_t GetIndex(ULong64_t mask)
 107   {
 108     for ( Int_t i = 0; i < 64; ++i )
 109     {
 110       if ( mask & ( 1ull << i ) ) return i+1;
 111     }
 112     return 0;
 113   }
 114
 115
 116   //______________________________________________________________________________
 117   Double_t Mu(Double_t L0B, Double_t Nb)
 118   {
 119     /// L0B = trigger rate before any veto
 120     /// Nb = number of crossing bunches
 121
 122     Double_t p0 = 1-L0B/(Nb*11245.0); // proba to get *no* collision per bunch crossing
 123
 124     return -TMath::Log(p0);
 125   }
 126
 127   //______________________________________________________________________________
 128   TCanvas* NewCanvas(const char* name)
 129   {
 130     TCanvas* c = new TCanvas(name,name);
 131     c->SetLeftMargin(0.15);
 132     return c;
 133   }
 134
 135   //______________________________________________________________________________
 136   void TimeAxis(TGraph* g)
 137   {
 138 //    g->GetXaxis()->SetTimeDisplay(1);
 139 //    g->GetXaxis()->SetTimeFormat("%d/%m %H:%M%F2010-12-31 24:00:00");
 140 //    g->GetXaxis()->SetNdivisions(505);
 141     g->GetYaxis()->SetTitleOffset(1.5);
 142
 143     g->GetXaxis()->SetTimeDisplay(1);
 144     //  g->GetXaxis()->SetTimeFormat("%d/%m %H:%M%F2010-12-31 24:00:00");
 145     g->GetXaxis()->SetTimeFormat("%d/%m %H:%M");
 146     g->GetXaxis()->SetTimeOffset(0,"gmt");
 147     g->GetXaxis()->SetNdivisions(505);
 148
 149   }
 150
 151 }
 152
 153 //_____________________________________________________________________________
 154 AliAnalysisTriggerScalers::AliAnalysisTriggerScalers(Int_t runNumber, const char* ocdbPath) : fRunList(), fOCDBPath(ocdbPath), fShouldCorrectForPileUp(kTRUE), fCrossSectionUnit("UB")
 155 {
 156   // ctor using a single run number
 157   SetRunList(runNumber);
 158 }
 159
 160 //_____________________________________________________________________________
 161 AliAnalysisTriggerScalers::AliAnalysisTriggerScalers(const std::vector<int>& runs, const char* ocdbPath) : fRunList(), fOCDBPath(ocdbPath), fShouldCorrectForPileUp(kTRUE), fCrossSectionUnit("UB")
 162 {
 163   // ctor using a vector of run numbers
 164   SetRunList(runs);
 165 }
 166
 167 //_____________________________________________________________________________
 168 AliAnalysisTriggerScalers::AliAnalysisTriggerScalers(const std::set<int>& runs, const char* ocdbPath) : fRunList(), fOCDBPath(ocdbPath), fShouldCorrectForPileUp(kTRUE), fCrossSectionUnit("UB")
 169 {
 170   // ctor using a set of run numbers
 171   SetRunList(runs);
 172 }
 173
 174 //_____________________________________________________________________________
 175 AliAnalysisTriggerScalers::AliAnalysisTriggerScalers(const char* runlist, const char* ocdbPath) : fRunList(), fOCDBPath(ocdbPath), fShouldCorrectForPileUp(kTRUE), fCrossSectionUnit("UB")
 176 {
 177   // ctor from a run list (txt file)
 178   SetRunList(runlist);
 179 }
 180
 181 //_____________________________________________________________________________
 182 AliAnalysisTriggerScalers::~AliAnalysisTriggerScalers()
 183 {
 184   // dtor
 185 }
 186
 187 //______________________________________________________________________________
 188 Bool_t AliAnalysisTriggerScalers::CheckRecord(const AliTriggerScalersRecord& record,
 189                                               UInt_t index,
 190                                               UInt_t refb,
 191                                               UInt_t refa,
 192                                               UInt_t timelapse) const
 193 {
 194   /// Check if a trigger scaler record should be kept for our
 195   /// luminosity computations
 196
 197   const AliTriggerScalers* scaler = record.GetTriggerScalersForClass(index);
 198
 199   UInt_t a = scaler->GetLOCA() - refa;
 200   UInt_t b = scaler->GetLOCB() - refb;
 201
 202   Bool_t ok(kTRUE);
 203
 204   if ( b <= 2 || ( a <= 2 ) || timelapse <= 9 ) // empirical cuts
 205   {
 206     ok = kFALSE;
 207   }
 208
 209   return ok;
 210 }
 211
 212
 213 //______________________________________________________________________________
 214 void AliAnalysisTriggerScalers::DrawFill(Int_t run1, Int_t run2, double ymin, double ymax, const char* label)
 215 {
 216   // Draw a yellow box to indicate a run range
 217
 218   AliDebugClass(1,Form("RUN1 %09d RUN2 %09d YMIN %e YMAX %e %s",
 219                        run1,run2,ymin,ymax,label));
 220   TBox* b = new TBox(run1*1.0,ymin,run2*1.0,ymax);
 221   b->SetFillColor(5);
 222   b->Draw();
 223   TText* text = new TText((run1+run2)/2.0,ymin + (ymax-ymin)*0.85,label);
 224   text->SetTextSize(0.025);
 225   text->SetTextFont(42);
 226   text->SetTextAlign(23);
 227   text->SetTextAngle(45);
 228   text->Draw();
 229 }
 230
 231 //_____________________________________________________________________________
 232 void AliAnalysisTriggerScalers::DrawFills(Double_t ymin, Double_t ymax)
 233 {
 234   /// Draw the fill ranges corresponding to the list of runs
 235   /// Note that this method will scan the OCDB to get the run -> fill number relationship,
 236   /// so it's better in this case to use a local copy of the OCDB if you have one. Otherwise
 237   /// it will be long.
 238   ///
 239   std::map<int, std::pair<int,int> > fills;
 240
 241   GetFillBoundaries(fills);
 242
 243   std::map<int, std::pair<int,int> >::const_iterator it;
 244
 245   for ( it = fills.begin(); it != fills.end(); ++it )
 246   {
 247     const std::pair<int,int>& p = it->second;
 248     TString fillnumber;
 249     fillnumber.Form("%d",it->first);
 250     DrawFill(p.first,p.second,ymin,ymax,fillnumber.Data());
 251   }
 252 }
 253
 254 //_____________________________________________________________________________
 255 void AliAnalysisTriggerScalers::GetCTPObjects(Int_t runNumber,
 256                                               AliTriggerConfiguration*& tc,
 257                                               AliTriggerRunScalers*& trs,
 258                                               AliLHCData*& lhc) const
 259 {
 260   /// For a given run, get the CTP objects we need to do our work
 261   ///
 262   /// FIXME: this is the method that should probably be virtualized so
 263   /// we can get those objects from either OCDB or run based container in AODs/ESDs
 264
 265   tc = static_cast<AliTriggerConfiguration*>(GetOCDBObject("GRP/CTP/Config",runNumber));
 266   trs = static_cast<AliTriggerRunScalers*>(GetOCDBObject("GRP/CTP/Scalers",runNumber));
 267   lhc = static_cast<AliLHCData*>(GetOCDBObject("GRP/GRP/LHCData",runNumber));
 268 }
 269
 270 //_____________________________________________________________________________
 271 void AliAnalysisTriggerScalers::GetFillBoundaries(std::map<int, std::pair<int,int> >& fills)
 272 {
 273   /// Given a list of runs get the fills they are in
 274
 275   fills.clear();
 276
 277   for ( std::vector<int>::const_iterator it = fRunList.begin(); it != fRunList.end(); ++it )
 278   {
 279     int runnumber = *it;
 280
 281     int fill = GetFillNumberFromRunNumber(runnumber);
 282
 283     if (fills.count(fill))
 284     {
 285       std::pair<int,int>& p = fills[fill];
 286       p.first = TMath::Min(runnumber,p.first);
 287       p.second = TMath::Max(runnumber,p.second);
 288     }
 289     else
 290     {
 291       fills[fill] = make_pair<int,int>(runnumber,runnumber);
 292     }
 293   }
 294 }
 295
 296
 297 //______________________________________________________________________________
 298 Int_t AliAnalysisTriggerScalers::GetFillNumberFromRunNumber(Int_t runNumber)
 299 {
 300   /// Get the fill number of a run
 301
 302   AliLHCData* lhcData = static_cast<AliLHCData*>(GetOCDBObject("GRP/GRP/LHCData",runNumber));
 303   if (lhcData)
 304   {
 305     Int_t fillNumber = lhcData->GetFillNumber();
 306     if ( fillNumber == 0 && runNumber == 189694)
 307     {
 308       // manual hack because GRP info incorrect for this run ?
 309       fillNumber = 3135;
 310     }
 311
 312     return fillNumber;
 313   }
 314   return -1;
 315 }
 316
 317 //______________________________________________________________________________
 318 TObject* AliAnalysisTriggerScalers::GetOCDBObject(const char* path, Int_t runNumber) const
 319 {
 320   // Helper method to get an object from the OCDB
 321
 322   if ( !AliCDBManager::Instance()->IsDefaultStorageSet() )
 323   {
 324     AliCDBManager::Instance()->SetDefaultStorage(fOCDBPath.Data());
 325   }
 326
 327   AliCDBManager::Instance()->SetRun(runNumber);
 328
 329   gErrorIgnoreLevel=kWarning; // to avoid all the TAlienFile::Open messages...
 330
 331   AliCDBEntry* e = AliCDBManager::Instance()->Get(path);
 332   return e ? e->GetObject() : 0x0;
 333 }
 334
 335
 336 //______________________________________________________________________________
 337 TString AliAnalysisTriggerScalers::GetLHCPeriodFromRunNumber(Int_t runNumber) const
 338 {
 339   /// Get a list of LHC periods from a list of run numbers
 340
 341   AliGRPObject* grp = static_cast<AliGRPObject*>(GetOCDBObject("GRP/GRP/Data",runNumber));
 342
 343   if (!grp) return "";
 344
 345   return grp->GetLHCPeriod();
 346 }
 347
 348 //_____________________________________________________________________________
 349 void
 350 AliAnalysisTriggerScalers::GetLuminosityTriggerAndCrossSection(Int_t runNumber,
 351                                                                TString& lumiTriggerClassName,
 352                                                                Double_t& lumiTriggerCrossSection,
 353                                                                Double_t& lumiTriggerCrossSectionError) const
 354 {
 355   /// For one given run, get the trigger class name to be used as the luminometer,
 356   /// and its corresponding cross-section
 357   ///
 358   /// FIXME: all is hard-coded for now... (use OADB for this ?)
 359
 360   lumiTriggerClassName="";
 361   lumiTriggerCrossSection=0.0;
 362   lumiTriggerCrossSectionError=0.0; // FIXME
 363
 364   TString period = GetLHCPeriodFromRunNumber(runNumber);
 365
 366   if ( period.BeginsWith("LHC11") )
 367   {
 368     AliError("Not implemented yet");
 369   }
 370   else if ( period.BeginsWith("LHC12") )
 371   {
 372      if ( period == "LHC12h" || period == "LHC12i" )
 373      {
 374        // pp 8 TeV main-satellites
 375        lumiTriggerClassName = "C0TVX-S-NOPF-ALLNOTRD";
 376        lumiTriggerCrossSection = 28.0; // FIXME: not from a vdM yet
 377      }
 378     else
 379     {
 380       AliError("Not implemented yet");
 381     }
 382   }
 383   else if ( period.BeginsWith("LHC13") )
 384   {
 385     if ( period == "LHC13g" )
 386     {
 387       // pp 2.76 TeV
 388       lumiTriggerClassName = "C0TVX-B-NOPF-ALLNOTRD";
 389       lumiTriggerCrossSection = 18.0; // FIXME: not from a vdM yet
 390     }
 391     else
 392     {
 393       // p-Pb or Pb-p 5.02 TeV
 394       lumiTriggerClassName = "C0TVX-B-NOPF-ALLNOTRD";
 395       lumiTriggerCrossSection = 0.755*2000; // FIXME: not from a vdM yet
 396     }
 397   }
 398   else
 399   {
 400     AliError("Not implemented yet");
 401   }
 402
 403   float csMult(1.0);
 404
 405   if ( fCrossSectionUnit=="PB" )
 406   {
 407     csMult=1E9;
 408   }
 409   else if (fCrossSectionUnit=="NB")
 410   {
 411     csMult=1E6;
 412   }
 413   else if ( fCrossSectionUnit=="UB" )
 414   {
 415     csMult=1E3;
 416   }
 417   else if ( fCrossSectionUnit=="MB" )
 418   {
 419     csMult=1.0;
 420   }
 421
 422   lumiTriggerCrossSection *= csMult;
 423   lumiTriggerCrossSectionError *= csMult;
 424 }
 425
 426 //_____________________________________________________________________________
 427 void AliAnalysisTriggerScalers::GetLHCPeriodBoundaries(std::map<std::string, std::pair<int,int> >& periods)
 428 {
 429   /// Given a list of runs get the fills they are in
 430
 431   periods.clear();
 432
 433   for ( std::vector<int>::const_iterator it = fRunList.begin(); it != fRunList.end(); ++it )
 434   {
 435     int runnumber = *it;
 436
 437     std::string period = GetLHCPeriodFromRunNumber(runnumber).Data();
 438
 439     if (periods.count(period))
 440     {
 441       std::pair<int,int>& p = periods[period];
 442       p.first = TMath::Min(runnumber,p.first);
 443       p.second = TMath::Max(runnumber,p.second);
 444     }
 445     else
 446     {
 447       periods[period] = make_pair<int,int>(runnumber,runnumber);
 448     }
 449   }
 450 }
 451
 452
 453 //______________________________________________________________________________
 454 Int_t AliAnalysisTriggerScalers::GetTriggerInput(Int_t runNumber, const char* inputname)
 455 {
 456   /// Get one trigger input for a given run
 457   AliTriggerConfiguration* tc = static_cast<AliTriggerConfiguration*>(GetOCDBObject("GRP/CTP/Config",runNumber));
 458   if (!tc) return -1;
 459
 460   const TObjArray& inputs = tc->GetInputs();
 461
 462   AliTriggerInput* ti = static_cast<AliTriggerInput*>(inputs.FindObject(inputname));
 463
 464   if (!ti) return -1;
 465
 466   return ti->GetSignature();
 467 }
 468
 469 //______________________________________________________________________________
 470 Float_t
 471 AliAnalysisTriggerScalers::GetPauseAndConfigCorrection(Int_t runNumber, const char* triggerClassName)
 472 {
 473   /// Tries to estimate the duration of the Pause and Configure phase(s) in a given run
 474   /// Probably not really accurate for the moment though.
 475
 476   ULong64_t total(0);
 477   ULong64_t nopac(0);
 478
 479   AliTriggerConfiguration* tc = static_cast<AliTriggerConfiguration*>(GetOCDBObject("GRP/CTP/Config",runNumber));
 480
 481   AliTriggerRunScalers* trs = static_cast<AliTriggerRunScalers*>(GetOCDBObject("GRP/CTP/Scalers",runNumber));
 482   const TObjArray* scalers = trs->GetScalersRecords();
 483
 484   TIter next(scalers);
 485   AliTriggerScalersRecord* record;
 486   UInt_t reft(0);
 487   UInt_t refl0b(0);
 488   UInt_t refl0a(0);
 489   Bool_t first(kTRUE);
 490
 491   Int_t classIndex = tc->GetClassIndexFromName(triggerClassName);
 492
 493   while ( ( record = static_cast<AliTriggerScalersRecord*>(next()) ) )
 494   {
 495     const AliTriggerScalers* scaler = record->GetTriggerScalersForClass(classIndex);
 496
 497     const AliTimeStamp* ats = record->GetTimeStamp();
 498
 499     UInt_t seconds = ats->GetSeconds();// - TTimeStamp::GetZoneOffset();
 500
 501     TTimeStamp ts(seconds,ats->GetMicroSecs());
 502
 503     UInt_t l0b = scaler->GetLOCB() - refl0b;
 504     UInt_t l0a = scaler->GetLOCA() - refl0a;
 505     UInt_t timelapse = seconds - reft;
 506
 507     if ( l0b <=2 || timelapse < 9 ) continue;
 508
 509     reft = seconds;
 510     refl0b = scaler->GetLOCB();
 511     refl0a = scaler->GetLOCA();
 512
 513     if ( first )
 514     {
 515       first = kFALSE;
 516       continue;
 517     }
 518
 519     total += l0b;
 520
 521     if ( l0a > 0 )
 522     {
 523       nopac += l0b;
 524     }
 525   }
 526
 527   return total > 0 ? nopac*1.0/total : 0.0;
 528 }
 529
 530 //______________________________________________________________________________
 531 AliAnalysisTriggerScalerItem*
 532 AliAnalysisTriggerScalers::GetTriggerScaler(Int_t runNumber, const char* level, const char* triggerClassName)
 533 {
 534   // Get the scaler for a given trigger class and a given trigger level.
 535   // Returned object must be deleted by the user.
 536
 537   AliTriggerConfiguration* tc = static_cast<AliTriggerConfiguration*>(GetOCDBObject("GRP/CTP/Config",runNumber));
 538   AliTriggerRunScalers* trs = static_cast<AliTriggerRunScalers*>(GetOCDBObject("GRP/CTP/Scalers",runNumber));
 539   AliGRPObject* grp = static_cast<AliGRPObject*>(GetOCDBObject("GRP/GRP/Data",runNumber));
 540
 541   if (!tc || !trs || !grp) return 0x0;
 542
 543   TString diCurrent(Form("L3:%5.0f;DIP:%5.0f [L3:%d;DIP:%d]",
 544                          grp->GetL3Current((AliGRPObject::Stats)0),
 545                          grp->GetDipoleCurrent((AliGRPObject::Stats)0),
 546                          grp->GetL3Polarity(),
 547                          grp->GetDipolePolarity()));
 548
 549   const TObjArray& trClasses = tc->GetClasses();
 550
 551   const TObjArray* scalers = trs->GetScalersRecords();
 552   const AliTriggerScalersRecord* begin = (AliTriggerScalersRecord*)(scalers->First());
 553   const AliTriggerScalersRecord* end = (AliTriggerScalersRecord*)(scalers->Last());
 554
 555   time_t duration = TMath::Nint((end->GetTimeStamp()->GetBunchCross() - begin->GetTimeStamp()->GetBunchCross())*AliTimeStamp::fNanosecPerBC*1E-9);
 556
 557   AliTriggerClass* triggerClass = static_cast<AliTriggerClass*>(trClasses.FindObject(triggerClassName));
 558   if (!triggerClass)
 559   {
 560     return 0x0;
 561   }
 562   UChar_t index = GetIndex(triggerClass->GetMask());
 563
 564   ULong64_t value(0);
 565
 566   const AliTriggerScalers* sbegin = begin->GetTriggerScalersForClass(index);
 567   const AliTriggerScalers* send = end->GetTriggerScalersForClass(index);
 568
 569   TString swhat(level);
 570   swhat.ToUpper();
 571
 572   if ( swhat.BeginsWith("L0A") )
 573   {
 574     value = send->GetLOCA() - sbegin->GetLOCA();
 575   }
 576   else if ( swhat.BeginsWith("L0B") )
 577   {
 578     value = send->GetLOCB() - sbegin->GetLOCB();
 579   }
 580   else if ( swhat.BeginsWith("L1A") )
 581   {
 582     value = send->GetL1CA() - sbegin->GetL1CA();
 583   }
 584   else if ( swhat.BeginsWith("L1B") )
 585   {
 586     value = send->GetL1CB() - sbegin->GetL1CB();
 587   }
 588   else if ( swhat.BeginsWith("L2A") )
 589   {
 590     value = send->GetL2CA() - sbegin->GetL2CA();
 591   }
 592   else if ( swhat.BeginsWith("L2B") )
 593   {
 594     value = send->GetL2CB() - sbegin->GetL2CB();
 595   }
 596   Int_t ds(1);
 597
 598   // FIXME: get the downscaling factor here for all cases (only BC1 supported here so far)
 599   if ( TString(triggerClassName).Contains("_B1") )
 600   {
 601     ds = GetTriggerInput(runNumber,"BC1");
 602     if (!ds) ds=1;
 603   }
 604
 605   swhat.ReplaceAll("RATE","");
 606
 607   return new AliAnalysisTriggerScalerItem(runNumber,swhat.Data(),diCurrent.Data(),triggerClass->GetName(),value,triggerClass->GetBCMask(),ds,duration);
 608 }
 609
 610 //______________________________________________________________________________
 611 TGraph*
 612 AliAnalysisTriggerScalers::IntegratedLuminosityGraph(const char* triggerClassName, const char* triggerClassNameForPACEstimate)
 613 {
 614   // Get the integrated luminosity graph for one trigger
 615
 616   if ( fRunList.empty() )
 617   {
 618     AliError("Cannot work without a run list");
 619     return 0x0;
 620   }
 621
 622   Double_t offset(0.0);
 623
 624   std::vector<double> vx;
 625   std::vector<double> vy;
 626
 627   Double_t x,y;
 628
 629   for ( std::vector<int>::size_type i = 0; i < fRunList.size(); ++i )
 630   {
 631     TGraph* g = IntegratedLuminosityGraph(fRunList[i],triggerClassName,triggerClassNameForPACEstimate);
 632
 633     if (!g) continue;
 634
 635     for ( Int_t j = 0 ; j < g->GetN(); ++j )
 636     {
 637       g->GetPoint(j,x,y);
 638       vx.push_back(x);
 639       vy.push_back(y+offset);
 640     }
 641
 642     offset += y;
 643
 644     delete g;
 645   }
 646
 647   TGraph* g = new TGraph(vx.size(),&vx[0],&vy[0]);
 648
 649   TimeAxis(g);
 650
 651   return g;
 652 }
 653
 654 //______________________________________________________________________________
 655 TGraph*
 656 AliAnalysisTriggerScalers::IntegratedLuminosityGraph(Int_t runNumber, const char* triggerClassName, const char* triggerClassNameForPACEstimate)
 657 {
 658   // Get the integrated luminosity graph for one trigger for the given run
 659   //
 660   // the triggerClassNameForPACEstimate is only used if triggerClassName is,
 661   // for the given run, the same as the luminometer trigger class name,
 662   // and should be the trigger class with the higher and non-zero L0A rate (except during PACs)
 663   //
 664
 665   AliTriggerConfiguration* tc(0x0);
 666   AliTriggerRunScalers* trs(0x0);
 667   AliLHCData* lhc(0x0);
 668
 669   GetCTPObjects(runNumber,tc,trs,lhc);
 670
 671   if (!tc || !trs || !lhc)
 672   {
 673     return 0x0;
 674   }
 675
 676   const TObjArray* scalerRecords= trs->GetScalersRecords();
 677   const TObjArray& triggerClasses = tc->GetClasses();
 678
 679   TString lumiTriggerClassName;
 680   Double_t lumiTriggerCrossSection(0.0);
 681   Double_t lumiTriggerCrossSectionError(0.0);
 682
 683   GetLuminosityTriggerAndCrossSection(runNumber,lumiTriggerClassName,lumiTriggerCrossSection,lumiTriggerCrossSectionError);
 684
 685   AliTriggerClass* lumiTriggerClass = static_cast<AliTriggerClass*>(triggerClasses.FindObject(lumiTriggerClassName));
 686
 687   if (!lumiTriggerClass)
 688   {
 689     AliError(Form("Could not find lumiTriggerClass=%s",lumiTriggerClassName.Data()));
 690     return 0x0;
 691   }
 692
 693   if (lumiTriggerCrossSection==0.0)
 694   {
 695     AliError(Form("Could not get cross-section for lumiTriggerClass=%s",lumiTriggerClassName.Data()));
 696     return 0x0;
 697   }
 698
 699   AliTriggerClass* triggerClass = static_cast<AliTriggerClass*>(triggerClasses.FindObject(triggerClassName));
 700
 701   if (!triggerClass)
 702   {
 703     AliError(Form("Could not find triggerClass=%s",triggerClassName));
 704     return 0x0;
 705   }
 706
 707   Int_t nbcx = NumberOfInteractingBunches(*lhc);
 708
 709   if (nbcx <= 0 && ShouldCorrectForPileUp())
 710   {
 711     AliError("Could not get the number of bunches, so won't be able to correct for pile-up");
 712   }
 713
 714   TIter nextScaler(scalerRecords);
 715   AliTriggerScalersRecord* record;
 716   UInt_t reft(0);
 717
 718   UInt_t refl0b[] = { 0, 0 };
 719   UInt_t refl0a[] = { 0, 0 };
 720
 721   UInt_t l0b[] = { 0, 0 };
 722   UInt_t l0a[] = { 0, 0 };
 723
 724   Bool_t first(kTRUE);
 725
 726   UChar_t classIndices[2];
 727
 728   // class 0 = class for luminomiter
 729   // class 1 = class for livetime estimate
 730   //           or for PAC estimate if triggerClassNameForPACEstimate
 731   //              is given and triggerClass is the same as the lumi class
 732
 733   classIndices[0] = GetIndex(lumiTriggerClass->GetMask());
 734   classIndices[1] = GetIndex(triggerClass->GetMask());
 735
 736   Bool_t sameClass = ( classIndices[0] == classIndices[1] );
 737   Bool_t pacRemoval(kFALSE);
 738
 739   if ( sameClass && strlen(triggerClassNameForPACEstimate) > 0 )
 740   {
 741     AliTriggerClass* triggerClassForPACEstimate = static_cast<AliTriggerClass*>(triggerClasses.FindObject(triggerClassNameForPACEstimate));
 742     if (!triggerClassForPACEstimate)
 743     {
 744       AliError(Form("Could not find triggerClassForPACEstimate=%s. Will not correct for PAC durations",triggerClassNameForPACEstimate));
 745       return 0x0;
 746     }
 747     classIndices[1] = GetIndex(triggerClassForPACEstimate->GetMask());
 748     sameClass = ( classIndices[0] == classIndices[1] );
 749     if (!sameClass)
 750     {
 751       pacRemoval=kTRUE;
 752     }
 753   }
 754
 755   ULong64_t counter(0);
 756
 757   std::vector<Double_t> vseconds;
 758   std::vector<Double_t> vlivetime;
 759   std::vector<Double_t> vlumi;
 760
 761   while ( ( record = static_cast<AliTriggerScalersRecord*>(nextScaler()) ) )
 762   {
 763     const AliTimeStamp* ats = record->GetTimeStamp();
 764
 765     UInt_t seconds = ats->GetSeconds();// - TTimeStamp::GetZoneOffset();
 766
 767     TTimeStamp ts(seconds,ats->GetMicroSecs());
 768
 769     UInt_t timelapse = seconds - reft;
 770
 771     Bool_t ok = sameClass || CheckRecord(*record,classIndices[1],refl0b[1],refl0a[1],timelapse);
 772
 773     if (ok) reft = seconds;
 774
 775     for ( Int_t i = 0; i < 2; ++i )
 776     {
 777       const AliTriggerScalers* scaler = record->GetTriggerScalersForClass(classIndices[i]);
 778
 779       if (ok)
 780       {
 781         l0b[i] = scaler->GetLOCB() - refl0b[i];
 782       }
 783       else
 784       {
 785         l0b[i] = 0;
 786       }
 787       l0a[i] = scaler->GetLOCA() - refl0a[i];
 788
 789       refl0b[i] = scaler->GetLOCB();
 790       refl0a[i] = scaler->GetLOCA();
 791     }
 792
 793     if ( first )
 794     {
 795       first = kFALSE;
 796       continue;
 797     }
 798
 799     Double_t factor(1.0);
 800
 801     if ( ShouldCorrectForPileUp() && nbcx > 0 && l0b[0] > 0 )
 802     {
 803       Double_t mu = Mu(l0b[0]/timelapse,nbcx);
 804       factor = mu/(1-TMath::Exp(-mu));
 805     }
 806
 807     vseconds.push_back(seconds);
 808
 809     Double_t lt(1.0);
 810
 811     if ( ok && !sameClass && !pacRemoval )
 812     {
 813       lt = (l0a[1]*1.0)/l0b[1];
 814     }
 815
 816     counter += TMath::Nint(factor*l0b[0]*lt);
 817
 818     vlumi.push_back( counter / lumiTriggerCrossSection );
 819
 820   }
 821
 822   if ( vseconds.empty() ) return 0x0;
 823
 824   TGraph* g = new TGraph(vseconds.size(),&vseconds[0],&vlumi[0]);
 825
 826   TimeAxis(g);
 827
 828   return g;
 829 }
 830
 831 //______________________________________________________________________________
 832 void AliAnalysisTriggerScalers::IntegratedLuminosity(const char* triggerList,
 833                                                      const char* lumiTrigger,
 834                                                      Double_t lumiCrossSection,
 835                                                      const char* csvOutputFile,
 836                                                      const char sep,
 837                                                      const char* csUnit)
 838 {
 839   // Compute the luminosity for a set of triggers
 840
 841   // for T0 based lumi (end of pp 2012), use lumiTrigger = C0TVX-S-NOPF-ALLNOTRD and lumiCrossSection = 28 mb (and csUnit="nb")
 842   //
 843   // for T0 based lumi (pPb 2013), use lumiTrigger = C0TVX-B-NOPF-ALLNOTRD
 844   // and lumiCrossSection = 0.755*2000 mb (and csUnit="mub")
 845
 846   // for T0 based lumi (pp 2.76 TeV 2013), use lumiTrigger = C0TVX-B-NOPF-ALLNOTRD
 847   // and lumiCrossSection = 18 mb (and csUnit="nb")
 848
 849   double intLumi(0);
 850
 851   std::map<std::string,float> lumiPerTrigger;
 852   std::map<int, std::map<std::string,float> > lumiPerFillPerTrigger;
 853
 854   std::map<std::string,time_t> durationPerTrigger;
 855
 856   TString sTriggerList(triggerList);
 857
 858   if ( sTriggerList.Length()==0 )
 859   {
 860     if ( lumiCrossSection < 30 && lumiCrossSection  > 20 )
 861     {
 862       // assuming it's pp 2012
 863       sTriggerList = "CMUL8-S-NOPF-ALLNOTRD,CMUL7-S-NOPF-ALLNOTRD,CMUL8-S-NOPF-MUON,CMUL7-S-NOPF-MUON";
 864
 865     // for C0MUL-SC-NOPF-MUON must use C0TVX-SC as lumiTrigger (with same cross-section as C0TVX=28mb)
 866     }
 867     else if ( lumiCrossSection < 30 )
 868     {
 869       // assuming it's pp 2013
 870       sTriggerList = "CMSL7-B-NOPF-MUON,CMSH7-B-NOPF-MUON,CMUL7-B-NOPF-MUON,CMSL7-B-NOPF-ALLNOTRD,CMSH7-B-NOPF-ALLNOTRD,CMUL7-B-NOPF-ALLNOTRD";
 871     }
 872     else
 873     {
 874       sTriggerList = "CMSL7-B-NOPF-MUON,CMSH7-B-NOPF-MUON,CMUL7-B-NOPF-MUON,CMSL7-B-NOPF-ALLNOTRD,CMSH7-B-NOPF-ALLNOTRD,CMUL7-B-NOPF-ALLNOTRD";
 875     }
 876   }
 877
 878   TObjArray* triggerArray = sTriggerList.Tokenize(",");
 879   triggerArray->SetOwner(kTRUE);
 880
 881   std::ofstream* out(0x0);
 882
 883   if ( TString(csvOutputFile).Length() > 0 )
 884   {
 885     out = new std::ofstream(gSystem->ExpandPathName(csvOutputFile));
 886     if (!out || out->bad())
 887     {
 888       delete out;
 889       out = 0x0;
 890     }
 891   }
 892
 893   TIter nextTrigger(triggerArray);
 894   TObjString* trigger;
 895
 896   if (out)
 897   {
 898     (*out) << "Fill number" << sep;
 899
 900     nextTrigger.Reset();
 901
 902     std::map<std::string,float>::const_iterator it;
 903
 904     while ( ( trigger = static_cast<TObjString*>(nextTrigger())))
 905     {
 906       lumiPerTrigger[trigger->String().Data()] = 0;
 907     }
 908
 909     for ( it = lumiPerTrigger.begin(); it != lumiPerTrigger.end(); ++it )
 910     {
 911       (*out) << "lumi from " << it->first.c_str() << sep;
 912     }
 913
 914     (*out) << "comments" << sep;
 915
 916     for ( it = lumiPerTrigger.begin(); it != lumiPerTrigger.end(); ++it )
 917     {
 918       (*out) << "recorded " <<  it->first.c_str() << " (integrated)" << sep;
 919     }
 920
 921     (*out) << Form("LHC delivered (%s^-1) per fill ",csUnit)
 922       << sep << Form("LHC delivered (%s^-1 integrated)",csUnit) << sep;
 923     (*out) << "lumi tot muon" << sep << "eff (%)" << sep;
 924     (*out) << std::endl;
 925
 926     nextTrigger.Reset();
 927   }
 928
 929   Int_t currentFillNumber(-1);
 930   Int_t fillNumber(0);
 931
 932   float csMult(1.0);
 933   TString scsUnit(csUnit);
 934   scsUnit.ToUpper();
 935
 936   if ( scsUnit=="PB" )
 937   {
 938     csMult=1E9;
 939   }
 940   else if (scsUnit=="NB")
 941   {
 942     csMult=1E6;
 943   }
 944   else if ( scsUnit=="UB" )
 945   {
 946     csMult=1E3;
 947   }
 948   else if ( scsUnit=="MB" )
 949   {
 950     csMult=1.0;
 951   }
 952   else
 953   {
 954     AliError(Form("Don't know how to deal with cross-section unit=%s",csUnit));
 955     return;
 956   }
 957
 958
 959   for ( std::vector<int>::size_type i = 0; i < fRunList.size(); ++i )
 960   {
 961     Int_t runNumber = fRunList[i];
 962     Bool_t atLeastOneTriggerFound(kFALSE);
 963
 964     // find out which trigger classes to use
 965
 966     AliTriggerConfiguration* tc = static_cast<AliTriggerConfiguration*>(GetOCDBObject("GRP/CTP/Config",runNumber));
 967     const TObjArray& trClasses = tc->GetClasses();
 968
 969     if (out)
 970     {
 971       fillNumber = GetFillNumberFromRunNumber(runNumber);
 972       if ( fillNumber == 0 )
 973       {
 974         AliError(Form("Got fillNumber = 0 for run %09d !",runNumber));
 975       }
 976
 977       if ( fillNumber != currentFillNumber )
 978       {
 979         std::map<std::string,float>::const_iterator it;
 980
 981         for ( it = lumiPerTrigger.begin(); it != lumiPerTrigger.end(); ++it )
 982         {
 983           lumiPerFillPerTrigger[fillNumber][it->first.c_str()] = 0;
 984         }
 985         currentFillNumber = fillNumber;
 986       }
 987     }
 988
 989     AliTriggerRunScalers* trs = static_cast<AliTriggerRunScalers*>(GetOCDBObject("GRP/CTP/Scalers",runNumber));
 990     const TObjArray* scalers = trs->GetScalersRecords();
 991     const AliTriggerScalersRecord* begin = (AliTriggerScalersRecord*)(scalers->First());
 992     const AliTriggerScalersRecord* end = (AliTriggerScalersRecord*)(scalers->Last());
 993
 994     time_t duration = TMath::Nint((end->GetTimeStamp()->GetBunchCross() - begin->GetTimeStamp()->GetBunchCross())*AliTimeStamp::fNanosecPerBC*1E-9);
 995
 996     nextTrigger.Reset();
 997
 998     TString lumiTriggerClassName(lumiTrigger);
 999     float lumiSigma = lumiCrossSection*csMult; //in csUnit
1000     AliAnalysisTriggerScalerItem* lumiB = GetTriggerScaler(runNumber,"L0B",lumiTriggerClassName.Data());
1001
1002     if (!lumiB)
1003     {
1004       AliError(Form("Did not find lumiTrigger %s for run %09d",lumiTriggerClassName.Data(),runNumber));
1005       continue;
1006     }
1007
1008     Float_t pacCorrection(1.0);
1009
1010     while ( ( trigger = static_cast<TObjString*>(nextTrigger()) ) )
1011     {
1012       TString muTriggerClassName(trigger->String());
1013       Int_t n(0);
1014
1015       if ( !trClasses.FindObject(muTriggerClassName.Data() ) )
1016       {
1017         continue;
1018       }
1019
1020       if ( muTriggerClassName.Contains("CMUL8") ) ++n;
1021       if ( muTriggerClassName.Contains("CMUL7") ) ++n;
1022
1023       if ( n>1 )
1024       {
1025         AliError(Form("More than 1 relevant trigger class found for run %09d ! Check that !",runNumber));
1026         trClasses.Print();
1027         continue;
1028       }
1029
1030       AliAnalysisTriggerScalerItem* muonA = GetTriggerScaler(runNumber,"L0A",muTriggerClassName.Data());
1031       AliAnalysisTriggerScalerItem* muonB = GetTriggerScaler(runNumber,"L0B",muTriggerClassName.Data());
1032
1033       if (!muonA || !muonB) continue;
1034
1035       atLeastOneTriggerFound = kTRUE;
1036
1037       Float_t ratio(0);
1038
1039       if ( muonB->ValueCorrectedForDownscale() > 0 )
1040       {
1041         ratio = muonA->ValueCorrectedForDownscale()/muonB->ValueCorrectedForDownscale();
1042       }
1043
1044       ratio *= lumiB->ValueCorrectedForDownscale()/lumiSigma;
1045
1046       if ( muTriggerClassName.BeginsWith("CMUL"))
1047       {
1048         intLumi += ratio;
1049       }
1050
1051       if ( muTriggerClassName.Contains("CMSH") )
1052       {
1053         pacCorrection = GetPauseAndConfigCorrection(runNumber,muTriggerClassName.Data());
1054       }
1055
1056       lumiPerTrigger[muTriggerClassName.Data()] += ratio;
1057       durationPerTrigger[muTriggerClassName.Data()] += duration;
1058       lumiPerFillPerTrigger[currentFillNumber][muTriggerClassName.Data()] += ratio;
1059
1060     }
1061
1062     lumiPerTrigger[lumiTriggerClassName.Data()] += lumiB->ValueCorrectedForDownscale()/lumiSigma;
1063     durationPerTrigger[lumiTriggerClassName.Data()] += duration;
1064     lumiPerFillPerTrigger[currentFillNumber][lumiTriggerClassName.Data()] += lumiB->ValueCorrectedForDownscale()/lumiSigma;
1065
1066     TString lumiPACCorrected(Form("%s(-PAC)",lumiTriggerClassName.Data()));
1067
1068     lumiPerTrigger[lumiPACCorrected.Data()] += pacCorrection*lumiB->ValueCorrectedForDownscale()/lumiSigma;
1069     durationPerTrigger[lumiPACCorrected.Data()] += duration;
1070     lumiPerFillPerTrigger[currentFillNumber][lumiPACCorrected.Data()] += pacCorrection*lumiB->ValueCorrectedForDownscale()/lumiSigma;
1071
1072     if (!atLeastOneTriggerFound && sTriggerList.Contains("CMUL") )
1073     {
1074       AliError(Form("Found no relevant trigger for run %09d",runNumber));
1075     }
1076   }
1077
1078   AliInfo(Form("Integrated lumi %7.4f %s^-1",intLumi,csUnit));
1079
1080   std::map<std::string,float>::const_iterator it;
1081
1082   for ( it = lumiPerTrigger.begin(); it != lumiPerTrigger.end(); ++it )
1083   {
1084     AliInfo(Form("Trigger %30s Lumi %10.4f %s^-1 duration %10ld s",it->first.c_str(),it->second,csUnit,durationPerTrigger[it->first]));
1085   }
1086
1087   if (out)
1088   {
1089     std::map<int, std::map<std::string, float> >::const_iterator fit;
1090
1091     lumiPerTrigger.clear();
1092
1093     for ( fit = lumiPerFillPerTrigger.begin(); fit != lumiPerFillPerTrigger.end(); ++fit )
1094     {
1095       int fill = fit->first;
1096       std::map<std::string,float>::const_iterator tit;
1097       (*out) << fill << sep;
1098
1099       for ( tit = fit->second.begin(); tit != fit->second.end(); ++tit )
1100       {
1101         (*out) << Form("%e",tit->second) << sep;
1102       }
1103
1104       (*out) << sep; // comment (empty)
1105
1106       for ( tit = fit->second.begin(); tit != fit->second.end(); ++tit )
1107       {
1108         lumiPerTrigger[tit->first] += tit->second;
1109
1110         (*out) <<  Form("%e",lumiPerTrigger[tit->first]) << sep;
1111       }
1112       (*out) << sep << "0" << sep << "0" << sep << "0" << sep << "0" << sep << std::endl; // LHC per fill, LHC integrated, lumi tot muon , efficiency
1113     }
1114   }
1115   //
1116   delete triggerArray;
1117 }
1118
1119 //______________________________________________________________________________
1120 TGraph* AliAnalysisTriggerScalers::MakeGraph(const std::vector<int>& vx,
1121                                              const std::vector<int>& vex,
1122                                              const std::vector<double>& vy,
1123                                              const std::vector<double>& vey)
1124 {
1125   /// Build a graph from a set of stl vectors
1126
1127   if ( ! ( vx.size() == vex.size() && vx.size() == vy.size() && vx.size() == vey.size() ) )
1128   {
1129     std::cerr << "incompatible sizes" << std::endl;
1130     return 0x0;
1131   }
1132
1133   double* x = new double[vx.size()];
1134   double* ex = new double[vx.size()];
1135   double* y = new double[vx.size()];
1136   double* ey = new double[vx.size()];
1137
1138   for ( size_t i = 0; i < vx.size(); ++i )
1139   {
1140     x[i] = vx[i];
1141     ex[i] = vex[i];
1142     y[i] = vy[i];
1143     ey[i] = vey[i];
1144   }
1145
1146   TGraph* g = new TGraphErrors(vx.size(),x,y,ex,ey);
1147
1148   TimeAxis(g);
1149
1150   delete[] x;
1151   delete[] y;
1152   delete[] ex;
1153   delete[] ey;
1154
1155   return g;
1156 }
1157
1158 //______________________________________________________________________________
1159 Int_t AliAnalysisTriggerScalers::NumberOfInteractingBunches(const AliLHCData& lhc) const
1160 {
1161   /// Extract the number of colliding bunches from the LHC data
1162
1163   Int_t numberOfInteractingBunches(0);
1164   Int_t numberOfInteractingBunchesMeasured(0);
1165   Int_t nIBM2(0);
1166
1167   int beam1(0);
1168   int beam2(1);
1169
1170   AliLHCDipValI* val = lhc.GetBunchConfigDeclared(beam1,0);
1171
1172   for ( Int_t i = 0; i < val->GetSizeTotal(); ++i )
1173   {
1174     if ( val->GetValue(i) < 0 ) ++numberOfInteractingBunches;
1175   }
1176
1177   AliLHCDipValI* valm = lhc.GetBunchConfigMeasured(beam1,0);
1178
1179   for ( Int_t i = 0; i < valm->GetSizeTotal(); ++i )
1180   {
1181     if ( valm->GetValue(i) < 0 ) ++numberOfInteractingBunchesMeasured;
1182   }
1183
1184   valm = lhc.GetBunchConfigMeasured(beam2,0);
1185
1186   for ( Int_t i = 0; i < valm->GetSizeTotal(); ++i )
1187   {
1188     if ( valm->GetValue(i) <= 0 ) ++nIBM2;
1189   }
1190
1191   if ( numberOfInteractingBunches != numberOfInteractingBunchesMeasured ||
1192        numberOfInteractingBunches != nIBM2 )
1193   {
1194     AliWarning(Form("Got some different number of interacting bunches here ! NumberOfInteractingBunches=%3d NumberOfInteractingBunchesMeasured=%3d NIBM2=%3d",
1195                   numberOfInteractingBunches,numberOfInteractingBunchesMeasured,nIBM2));
1196   }
1197
1198   return numberOfInteractingBunches;
1199 }
1200
1201
1202 //______________________________________________________________________________
1203 TGraph* AliAnalysisTriggerScalers::PlotTrigger(const char* triggerClassName,
1204                                                const char* what)
1205 {
1206   // Plot one of the scalers (L0A,L0B,L0AOVERB,etc...) of a given triggerClass
1207   // Get one value per run.
1208
1209   std::vector<Float_t> x;
1210   std::vector<Float_t> y;
1211
1212   double integral(0);
1213   double mean(0);
1214
1215   for ( std::vector<int>::size_type i = 0; i < fRunList.size(); ++i )
1216   {
1217     Int_t runNumber = fRunList[i];
1218
1219     AliAnalysisTriggerScalerItem* s = GetTriggerScaler(runNumber,what,triggerClassName);
1220
1221     if (!s) continue;
1222
1223     x.push_back(runNumber);
1224
1225     Double_t value = s->ValueCorrectedForDownscale();
1226
1227     if ( TString(what).Contains("RATE") )
1228     {
1229       value = s->Rate();
1230     }
1231
1232     integral += value;
1233     mean += value;
1234
1235     y.push_back(value);
1236   }
1237
1238   if ( fRunList.size() ) {
1239     mean /= fRunList.size();
1240   }
1241
1242   AliInfo(Form("Integral %e mean %e",integral,mean));
1243
1244   return new TGraph(x.size(),&x[0],&y[0]);
1245 }
1246
1247 //______________________________________________________________________________
1248 TGraph* AliAnalysisTriggerScalers::PlotTriggerEvolution(const char* triggerClassName,
1249                                                         const char* what,
1250                                                         bool draw,
1251                                                         double* mean,
1252                                                         bool removeZeros)
1253 {
1254   /// Make a graph of "what" for a given trigger class.
1255   ///
1256   /// What can be :
1257   ///
1258   /// - L0B (level 0 before any veto applied)
1259   /// - L0A (level 0 after all vetoes)
1260   /// - L0AOVERB (L0A/L0B)
1261   /// - mu ( = -TMath::Log( 1 - P(0) ) where P(0) is the proba to have zero collisions in a bunch crossing)
1262   /// - pileupfactor = mu/(1-exp(-mu)) : the factor to apply to correct the cint1b count rate
1263   /// - vsnb = L0B/(NumberOfInteractingBunches*11245)
1264   /// - NBCX = NumberOfInteractingBunches
1265
1266   TString swhat(what);
1267   swhat.ToUpper();
1268
1269   if ( swhat.Contains(";"))
1270   {
1271     swhat.ReplaceAll(";",",");
1272     AliWarningClass("; is not a valid separator, replaced it with ,");
1273   }
1274
1275   int color(1);
1276   int marker(20);
1277
1278   TObjArray* a = swhat.Tokenize(",");
1279   if (a->GetEntries()>1)
1280   {
1281     TObjArray graphs;
1282     TIter next(a);
1283     TObjString* str(0x0);
1284     Double_t ymin(TMath::Limits<Double_t>::Max());
1285     Double_t ymax(0);
1286     Double_t xmin(TMath::Limits<Double_t>::Max());
1287     Double_t xmax(0);
1288     TGraph* g(0x0);
1289
1290     while ( ( str = static_cast<TObjString*>(next())) )
1291     {
1292       g = PlotTriggerEvolution(triggerClassName,str->String().Data(),false,0x0,removeZeros);
1293       graphs.Add(g);
1294       for ( Int_t i = 0; i < g->GetN(); ++i )
1295       {
1296         ymin = TMath::Min(ymin,g->GetY()[i]);
1297         ymax = TMath::Max(ymax,g->GetY()[i]);
1298       }
1299       xmin = TMath::Min(xmin,g->GetX()[0]);
1300       xmax = TMath::Max(xmax,g->GetX()[g->GetN()-1]);
1301     }
1302
1303     gStyle->SetOptTitle(0);
1304
1305     AliInfoClass(Form("x %e ; %e y %e ; %e",xmin,xmax,ymin,ymax));
1306     TH2* h = new TH2F("h",triggerClassName,100,xmin,xmax,100,ymin,ymax);
1307     h->SetStats(kFALSE);
1308     h->GetXaxis()->SetTimeDisplay(1);
1309     h->GetXaxis()->SetTimeFormat("%d/%m %H:%M");
1310     h->GetXaxis()->SetTimeOffset(0,"gmt");
1311     h->GetXaxis()->SetNdivisions(505);
1312     h->Draw();
1313
1314     TIter nextGraph(&graphs);
1315
1316     while ( ( g = static_cast<TGraph*>(nextGraph())) )
1317     {
1318       AliInfoClass(g->GetTitle());
1319       g->Draw("lp");
1320       g->SetLineColor(color);
1321       g->SetMarkerColor(color);
1322       g->SetMarkerStyle(marker);
1323       ++color;
1324       ++marker;
1325     }
1326     delete a;
1327     return 0x0;
1328   }
1329   delete a;
1330
1331   std::vector<int> vx;
1332   std::vector<int> vex;
1333   std::vector<double> vy;
1334   std::vector<double> vey;
1335
1336
1337   if (mean) *mean=0;
1338   double nvalues(0);
1339
1340   for ( std::vector<int>::size_type iRun = 0; iRun < fRunList.size(); ++iRun )
1341   {
1342     Int_t runNumber = fRunList[iRun];
1343
1344     AliTriggerConfiguration* tc(0x0);
1345     AliTriggerRunScalers* trs(0x0);
1346     AliLHCData* lhc(0x0);
1347
1348     GetCTPObjects(runNumber,tc,trs,lhc);
1349
1350     Int_t numberOfInteractingBunches = NumberOfInteractingBunches(*lhc);
1351
1352     const TObjArray* scalers = trs->GetScalersRecords();
1353
1354     const TObjArray& trClasses = tc->GetClasses();
1355     TIter next(&trClasses);
1356     AliTriggerClass* triggerClass;
1357
1358     while ( ( triggerClass = static_cast<AliTriggerClass*>(next()) ) )
1359     {
1360       UChar_t index = GetIndex(triggerClass->GetMask());
1361
1362       if ( !TString(triggerClass->GetName()).Contains(triggerClassName) ) continue;
1363
1364       TIter nextScaler(scalers);
1365       AliTriggerScalersRecord* record;
1366       UInt_t reft(0);
1367       UInt_t refl0b(0);
1368       UInt_t refl1b(0);
1369       UInt_t refl2b(0);
1370       UInt_t refl0a(0);
1371       UInt_t refl1a(0);
1372       UInt_t refl2a(0);
1373       Bool_t first(kTRUE);
1374
1375       while ( ( record = static_cast<AliTriggerScalersRecord*>(nextScaler()) ) )
1376       {
1377         const AliTriggerScalers* scaler = record->GetTriggerScalersForClass(index);
1378
1379
1380         const AliTimeStamp* ats = record->GetTimeStamp();
1381
1382         UInt_t seconds = ats->GetSeconds();// - TTimeStamp::GetZoneOffset();
1383
1384         TTimeStamp ts(seconds,ats->GetMicroSecs());
1385
1386         UInt_t l0b = scaler->GetLOCB() - refl0b;
1387         UInt_t l0a = scaler->GetLOCA() - refl0a;
1388         UInt_t l1b = scaler->GetL1CB() - refl1b;
1389         UInt_t l1a = scaler->GetL1CA() - refl1a;
1390         UInt_t l2b = scaler->GetL2CB() - refl2b;
1391         UInt_t l2a = scaler->GetL2CA() - refl2a;
1392         UInt_t timelapse = seconds - reft;
1393
1394         if ( removeZeros && ( l0b <= 2 || ( l0a <= 2 && l0a != 0 ) || timelapse <= 9 ) )
1395         {
1396           AliInfo(Form("Skipping point for %s l0b %d l0a %d timelapse %d ts %s",
1397                        triggerClassName,l0b,l0a,timelapse,ts.AsString()));
1398           continue;
1399         }
1400
1401         reft = seconds;
1402         refl0b = scaler->GetLOCB();
1403         refl0a = scaler->GetLOCA();
1404         refl1b = scaler->GetL1CB();
1405         refl1a = scaler->GetL1CA();
1406         refl2b = scaler->GetL2CB();
1407         refl2a = scaler->GetL2CA();
1408
1409         if ( first )
1410         {
1411           first = kFALSE;
1412           continue;
1413         }
1414
1415         double value(1.0);
1416         double error(0.0);
1417
1418         if (swhat.Contains("L0AOVERB") )
1419         {
1420           value = l0a*1.0/l0b;
1421           if ( l0a > 0 )
1422           {
1423             error = value*TMath::Sqrt(1.0/l0b+1.0/l0a);
1424           }
1425           else
1426           {
1427             error = 0.0;
1428           }
1429         }
1430         else if ( swhat.Contains("L0B") )
1431         {
1432           value = l0b;
1433         }
1434         else if (swhat.Contains("L0A") )
1435         {
1436           value = l0a;
1437         }
1438         else if ( swhat.Contains("L1B") )
1439         {
1440           value = l1b;
1441         }
1442         else if (swhat.Contains("L1A") )
1443         {
1444           value = l1a;
1445         }
1446         else if ( swhat.Contains("L2B") )
1447         {
1448           value = l2b;
1449         }
1450         else if (swhat.Contains("L2A") )
1451         {
1452           value = l2a;
1453         }
1454         else if ( swhat.Contains("MU") )
1455         {
1456           value = Mu(l0b/timelapse,numberOfInteractingBunches);
1457           error = 0.0; // FIXME
1458         }
1459         else if ( swhat.Contains("PILEUPFACTOR") )
1460         {
1461           Double_t mu = Mu(l0b/timelapse,numberOfInteractingBunches);
1462           value = mu/(1-TMath::Exp(-mu));
1463           error = -1.0; // FIXME
1464         }
1465         else if ( swhat.Contains("VSNB") )
1466         {
1467           value = l0b/(11245.0*numberOfInteractingBunches);
1468           error = -1.0; // FIXME
1469         }
1470         else if ( swhat.Contains("NBCX"))
1471         {
1472           value = numberOfInteractingBunches;
1473           error = 1.0; // FIXME
1474         }
1475         else
1476         {
1477           value = timelapse;
1478           AliInfo(Form("Unknown what %s",what));
1479         }
1480
1481         if ( ! swhat.Contains("OVER") && ! swhat.Contains("RATIO") &&
1482             ! swhat.Contains("MU") && ! swhat.Contains("PILEUPFACTOR") &&
1483             ! swhat.Contains("RAW") & ! swhat.Contains("NBCX") )
1484         {
1485           value /= timelapse;
1486         }
1487
1488         if ( !TMath::Finite(value) ) continue;
1489
1490         if ( !swhat.Contains("L0AOVERB") && error > 0 )
1491         {
1492           error = value >0 ? 1.0/TMath::Sqrt(value) : 0.0;
1493         }
1494
1495         if (removeZeros && TMath::Abs(value) < 1E-6 )
1496         {
1497           continue;
1498         }
1499         if (mean)
1500         {
1501           (*mean) += value;
1502           nvalues += 1.0;
1503         }
1504
1505         vx.push_back(seconds);
1506         vex.push_back(1);
1507
1508         vy.push_back(value);
1509
1510         vey.push_back(error);
1511
1512       }
1513
1514     }
1515
1516   }
1517
1518   if (mean && nvalues)
1519   {
1520     (*mean) /= nvalues;
1521   }
1522
1523   if ( vx.empty() ) return 0;
1524
1525   TGraph* g = MakeGraph(vx,vex,vy,vey);
1526   TString title(Form("TriggerEvolution%s-%s",triggerClassName,swhat.Data()));
1527
1528   g->SetName(title.Data());
1529   g->SetTitle(title.Data());
1530
1531   g->GetYaxis()->SetTitle(title.Data());
1532
1533   if (draw)
1534   {
1535     TCanvas* c = NewCanvas(g->GetName());
1536     g->SetLineColor(color);
1537     g->SetMarkerColor(color);
1538     g->SetMarkerStyle(marker);
1539     g->Draw("ALP");
1540     c->SaveAs(Form("%s.png",title.Data()));
1541   }
1542
1543
1544   return g;
1545 }
1546
1547 //______________________________________________________________________________
1548 TGraph* AliAnalysisTriggerScalers::PlotTriggerRatio(const char* triggerClassName1,
1549                                                     const char* what1,
1550                                                     const char* triggerClassName2,
1551                                                     const char* what2)
1552 {
1553   // Plot the ratio of two scalers.
1554   // Get one value per run.
1555
1556   std::vector<Float_t> x;
1557   std::vector<Float_t> y;
1558
1559   for ( std::vector<int>::size_type i = 0; i < fRunList.size(); ++i )
1560   {
1561     Int_t runNumber = fRunList[i];
1562
1563     AliAnalysisTriggerScalerItem* s1 = GetTriggerScaler(runNumber,what1,triggerClassName1);
1564     AliAnalysisTriggerScalerItem* s2 = GetTriggerScaler(runNumber,what2,triggerClassName2);
1565
1566     if (!s1 || !s2) continue;
1567
1568     x.push_back(runNumber);
1569
1570     Float_t ratio(0);
1571
1572     if ( s2->ValueCorrectedForDownscale() > 0 )
1573     {
1574       ratio = s1->ValueCorrectedForDownscale()/s2->ValueCorrectedForDownscale();
1575     }
1576
1577     y.push_back(ratio);
1578
1579       AliDebug(1,Form("RUN %09d %20s (%s) %12llu (%5d) %20s (%s) %12llu (%5d) R %7.2f",
1580                    runNumber,
1581                    triggerClassName1,what1,
1582                    s1->Value(),s1->DownscalingFactor(),
1583                    triggerClassName2,what2,
1584                    s2->Value(),s2->DownscalingFactor(),
1585                    ratio));
1586   }
1587
1588   return new TGraph(x.size(),&x[0],&y[0]);
1589 }
1590
1591 //______________________________________________________________________________
1592 TGraph* AliAnalysisTriggerScalers::PlotTriggerRatioEvolution(const char* triggerClassName1,
1593                                                              const char* what1,
1594                                                              const char* triggerClassName2,
1595                                                              const char* what2)
1596 {
1597   /// Plots the evolution of one trigger ratio
1598
1599   Bool_t draw(kFALSE);
1600   Bool_t removeZeros(kFALSE);
1601
1602   TGraph* g1 = PlotTriggerEvolution(triggerClassName1,what1,draw,0x0,removeZeros);
1603   TGraph* g2 = PlotTriggerEvolution(triggerClassName2,what2,draw,0x0,removeZeros);
1604
1605   if (!g1 || !g2) return 0x0;
1606
1607   if ( g1->GetN() != g2->GetN() )
1608   {
1609     AliError("Oups. Did not get the same number of points for the 2 graphs ?!");
1610     return 0x0;
1611   }
1612
1613   Double_t x1,x2;
1614   Double_t y1,y2;
1615   Double_t x1err,x2err;
1616   Double_t y1err,y2err;
1617
1618   TGraphErrors* g = new TGraphErrors(g1->GetN());
1619   Int_t j(0);
1620
1621   for ( Int_t i = 0; i < g1->GetN(); ++i )
1622   {
1623     g1->GetPoint(i,x1,y1);
1624     g2->GetPoint(i,x2,y2);
1625
1626     x1err = g1->GetErrorX(i);
1627     x2err = g2->GetErrorX(i);
1628
1629     y1err = g1->GetErrorY(i);
1630     y2err = g2->GetErrorY(i);
1631
1632     if  (x1!=x2)
1633     {
1634       AliError(Form("Points at %d don't have the same x ! ",i));
1635       continue;
1636     }
1637
1638     Double_t y(0.0);
1639
1640     if ( TMath::Abs(y2) > 1E-12 )
1641     {
1642       y = y1/y2;
1643     }
1644
1645     Double_t yerr(0.0);
1646
1647     if ( TMath::Abs(x1) > 1E-12)
1648     {
1649       yerr += TMath::Sqrt(x1err*x1err/(x1*x1));
1650     }
1651
1652     if ( TMath::Abs(x2) > 1E-12)
1653     {
1654       yerr += TMath::Sqrt(x2err*x2err/(x2*x2));
1655     }
1656
1657     yerr *= y;
1658
1659     g->SetPoint(j,x1,y);
1660     g->SetPointError(j,x1err,yerr);
1661
1662     ++j;
1663   }
1664
1665   delete g1;
1666   delete g2;
1667
1668   TimeAxis(g);
1669
1670   return g;
1671 }
1672
1673 //______________________________________________________________________________
1674 void AliAnalysisTriggerScalers::Print(Option_t* /* opt */) const
1675 {
1676   /// print our runlist
1677   AliAnalysisTriggerScalers::PrintIntegers(fRunList,',');
1678 }
1679
1680 //______________________________________________________________________________
1681 void AliAnalysisTriggerScalers::PrintIntegers(const std::vector<int>& integers,
1682                                               const char sep,
1683                                               std::ostream& out)
1684 {
1685   /// print a list of integers
1686   for ( std::vector<int>::size_type i = 0; i < integers.size(); ++i )
1687   {
1688     out << integers[i] << sep;
1689   }
1690   out << std::endl;
1691 }
1692
1693 //______________________________________________________________________________
1694 void AliAnalysisTriggerScalers::ReadIntegers(const char* filename,
1695                                              std::vector<int>& integers,
1696                                              Bool_t resetVector)
1697 {
1698   /// Read integers from filename, where integers are either
1699   /// separated by "," or by return carriage
1700     std::ifstream in(gSystem->ExpandPathName(filename));
1701   int i;
1702
1703   std::set<int> runset;
1704
1705   if (!resetVector)
1706   {
1707     for ( std::vector<int>::size_type j = 0; j < integers.size(); ++ j )
1708     {
1709       runset.insert(integers[j]);
1710     }
1711   }
1712
1713   char line[10000];
1714
1715   in.getline(line,10000,'\n');
1716
1717   TString sline(line);
1718
1719   if (sline.Contains(","))
1720   {
1721     TObjArray* a = sline.Tokenize(",");
1722     TIter next(a);
1723     TObjString* s;
1724     while ( ( s = static_cast<TObjString*>(next()) ) )
1725     {
1726       runset.insert(s->String().Atoi());
1727     }
1728     delete a;
1729   }
1730   else
1731   {
1732     runset.insert(sline.Atoi());
1733
1734     while ( in >> i )
1735     {
1736       runset.insert(i);
1737     }
1738   }
1739
1740   integers.clear();
1741
1742   for ( std::set<int>::const_iterator it = runset.begin(); it != runset.end(); ++it )
1743   {
1744     integers.push_back((*it));
1745   }
1746
1747   std::sort(integers.begin(),integers.end());
1748 }
1749
1750
1751 //______________________________________________________________________________
1752 void AliAnalysisTriggerScalers::SetRunList(Int_t runNumber)
1753 {
1754   // Make the runlist be a single run
1755   fRunList.clear();
1756   fRunList.push_back(runNumber);
1757 }
1758
1759 //______________________________________________________________________________
1760 void AliAnalysisTriggerScalers::SetRunList(const std::vector<int>& runs)
1761 {
1762   // Make the runlist be a single run
1763   fRunList = runs;
1764 }
1765
1766 //______________________________________________________________________________
1767 void AliAnalysisTriggerScalers::SetRunList(const std::set<int>& runs)
1768 {
1769   // Make the runlist be a single run
1770   fRunList.clear();
1771   for ( std::set<int>::const_iterator it = runs.begin(); it != runs.end(); ++it )
1772   {
1773     fRunList.push_back(*it);
1774   }
1775 }
1776
1777 //______________________________________________________________________________
1778 void AliAnalysisTriggerScalers::SetRunList(const char* runlist)
1779 {
1780   // Read the runlist from an ASCII file or a comma separated list
1781   // or a space separated list
1782
1783   fRunList.clear();
1784
1785   if ( TString(runlist).Contains(",") || TString(runlist).Contains(" ") )
1786   {
1787     TObjArray* runs = 0x0;
1788     if ( TString(runlist).Contains(",") )
1789     {
1790       runs = TString(runlist).Tokenize(",");
1791     }
1792     else
1793     {
1794       runs = TString(runlist).Tokenize(" ");
1795     }
1796     TIter next(runs);
1797     TObjString* s;
1798     std::set<int> runset;
1799
1800     while ( ( s = static_cast<TObjString*>(next()) ) )
1801     {
1802       runset.insert(s->String().Atoi());
1803     }
1804
1805     for ( std::set<int>::const_iterator it = runset.begin(); it != runset.end(); ++it )
1806     {
1807       fRunList.push_back((*it));
1808     }
1809
1810     std::sort(fRunList.begin(),fRunList.end());
1811
1812     delete runs;
1813   }
1814   else
1815   {
1816     ReadIntegers(runlist,fRunList);
1817   }
1818 }
1819
1820 //______________________________________________________________________________
1821 //______________________________________________________________________________
1822 //______________________________________________________________________________
1823 //______________________________________________________________________________
1824 //______________________________________________________________________________
1825 //______________________________________________________________________________
1826 //______________________________________________________________________________
1827
1828 ClassImp(AliAnalysisTriggerScalerItem)
1829
1830 //______________________________________________________________________________
1831 const char*
1832 AliAnalysisTriggerScalerItem::BCMaskName() const
1833 {
1834   // return bc mask name
1835   if ( BCMask() ) return BCMask()->GetName(); else return "";
1836 }
1837
1838 //______________________________________________________________________________
1839 Int_t AliAnalysisTriggerScalerItem::Compare(const TObject* obj) const
1840 {
1841   /// compare two scaler items (by means of their run number only)
1842   const AliAnalysisTriggerScalerItem* s = static_cast<const AliAnalysisTriggerScalerItem*>(obj);
1843
1844   if ( s->RunNumber() < RunNumber() )
1845   {
1846     return -1;
1847   }
1848   else if ( s->RunNumber() > RunNumber() )
1849   {
1850     return 1;
1851   }
1852   return 0;
1853 }
1854
1855 //______________________________________________________________________________
1856 void AliAnalysisTriggerScalerItem::Print(Option_t* opt) const
1857 {
1858   /// Printout
1859   TString sopt(opt);
1860
1861   if ( RunNumber() > 0 )
1862   {
1863     sopt.ToUpper();
1864
1865     if ( sopt.Contains("FULL") )
1866     {
1867     }
1868     else
1869     {
1870       std::cout << Form("RUN %6d CLASS %24s (%5s %4d BCs) SCALER %s %12llu DS %6d DURATION %ld",
1871                    RunNumber(),TriggerClassName(),
1872                    BCMaskName(),
1873                    BCMask() ? BCMask()->GetNUnmaskedBCs() : 0,
1874                    Level(),
1875                    Value(),DownscalingFactor(),Duration()) << std::endl;
1876     }
1877   }
1878   else
1879   {
1880     std::cout << Form("CLASS %24s SCALER %20llu NRUNS %d",
1881                  TriggerClassName(),Value(),NofRuns()) << std::endl;
1882   }
1883 }