ANALYSIS/AliPhysicsSelection.cxx

   1 /* $Id: AliPhysicsSelection.cxx 35782 2009-10-22 11:54:31Z jgrosseo $ */
   2
   3 /**************************************************************************
   4  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   5  *                                                                        *
   6  * Author: The ALICE Off-line Project.                                    *
   7  * Contributors are mentioned in the code where appropriate.              *
   8  *                                                                        *
   9  * Permission to use, copy, modify and distribute this software and its   *
  10  * documentation strictly for non-commercial purposes is hereby granted   *
  11  * without fee, provided that the above copyright notice appears in all   *
  12  * copies and that both the copyright notice and this permission notice   *
  13  * appear in the supporting documentation. The authors make no claims     *
  14  * about the suitability of this software for any purpose. It is          *
  15  * provided "as is" without express or implied warranty.                  *
  16  **************************************************************************/
  17
  18 //-------------------------------------------------------------------------
  19 //                      Implementation of   Class AliPhysicsSelection
  20 // This class selects collision candidates from data runs, applying selection cuts on triggers
  21 // and background rejection based on the content of the ESD
  22 //
  23 // Usage:
  24 //
  25 // Create the object:
  26 //   fPhysicsSelection = new AliPhysicsSelection;
  27 //
  28 // For MC data, call
  29 //   fPhysicsSelection->SetAnalyzeMC()
  30 //
  31 // To check if an event is a collision candidate, use:
  32 //   fPhysicsSelection->IsCollisionCandidate(fESD)
  33 //
  34 // After processing save the resulting histograms to a file with (a folder physics_selection
  35 //   will be created that contains the histograms):
  36 //   fPhysicsSelection->SaveHistograms("physics_selection")
  37 //
  38 // To print statistics after processing use:
  39 //   fPhysicsSelection->Print();
  40 //
  41 // The BX ids corresponding to real bunches crossings p2 are
  42 // automatically selected. You cannot process runs with different
  43 // filling schemes if this option is set. If you want to disable this,
  44 // use:
  45 //   fPhysicsSelection->SetUseBXNumbers(0);
  46 //
  47 //
  48 // If you are analizing muons and you want to keep the muon triggers
  49 // besides the CINT1B you can set:
  50 //   fPhysicsSelection->SetUseMuonTriggers();
  51 //
  52 //
  53 // To compute the Background automatically using the control triggers
  54 // use:
  55 //   fPhysicsSelection->SetComputeBG();
  56 // this will show the value of the Beam Gas, accidentals and good
  57 // events as additional rows in the statistic tables, but it will NOT
  58 // subtract the background automatically.
  59 // This option will only work for runs taken with the CINT1
  60 // suite. This options enables automatically also the usage of BX
  61 // numbers. You can only process one run at a time if you require this
  62 // options, because it uses the bunch intensity estimated run by run.
  63 //
  64 // The BG will usually be more important in the so-called "bin 0": the
  65 // class can also compute the statistics table for events in this
  66 // bin. Since the definition of bin 0 may in general change from
  67 // analysis to analysis, the user needs to provide a callback
  68 // implementing the definition of bin zero. The callback should be
  69 // implemented as a method in the analysis task and should override
  70 // the IsEventInBinZero method of AliAnalysisTaskSE, and should thus
  71 // have the the following prototype:
  72 //   Bool_t IsEventInBinZero();
  73 // It should return true if the event is in the bin 0 and it is set by
  74 // passing to the physics selection the NAME of the task where the
  75 // callback is implemented:
  76 //   fPhysicsSelection->SetBin0Callback("MyTask").
  77 //
  78 //
  79 // Usually the class selects the trigger scheme by itself depending on the run number.
  80 // Nevertheless, you can do that manually by calling AddCollisionTriggerClass() and AddBGTriggerClass()
  81 // Example:
  82 // To define the class CINT1B-ABCE-NOPF-ALL as collision trigger (those will be accepted as
  83 // collision candidates when they pass the selection):
  84 //   AddCollisionTriggerClass("+CINT1B-ABCE-NOPF-ALL #769 #3119");
  85 // To select on bunch crossing IDs in addition, use:
  86 //   AddCollisionTriggerClass("+CINT1B-ABCE-NOPF-ALL #769 #3119");
  87 // To define the class CINT1A-ABCE-NOPF-ALL as a background trigger (those will only be counted
  88 // for the control histograms):
  89 //   AddBGTriggerClass("+CINT1A-ABCE-NOPF-ALL");
  90 // You can also specify more than one trigger class in a string or you can require that some are *not*
  91 // present. The following line would require CSMBA-ABCE-NOPF-ALL, but CSMBB-ABCE-NOPF-ALL is not allowed
  92 // to be present:
  93 //   AddBGTriggerClass("+CSMBA-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL");
  94 //
  95 //   Origin: Jan Fiete Grosse-Oetringhaus, CERN
  96 //           Michele Floris, CERN
  97 //-------------------------------------------------------------------------
  98
  99 #include <Riostream.h>
 100 #include <TH1F.h>
 101 #include <TH2F.h>
 102 #include <TList.h>
 103 #include <TIterator.h>
 104 #include <TDirectory.h>
 105 #include <TObjArray.h>
 106
 107 #include <AliPhysicsSelection.h>
 108
 109 #include <AliTriggerAnalysis.h>
 110 #include <AliLog.h>
 111
 112 #include <AliESDEvent.h>
 113 #include <AliAnalysisTaskSE.h>
 114 #include "AliAnalysisManager.h"
 115 #include "TPRegexp.h"
 116
 117 ClassImp(AliPhysicsSelection)
 118
 119 AliPhysicsSelection::AliPhysicsSelection() :
 120   AliAnalysisCuts("AliPhysicsSelection", "AliPhysicsSelection"),
 121   fCurrentRun(-1),
 122   fMC(kFALSE),
 123   fCollTrigClasses(),
 124   fBGTrigClasses(),
 125   fTriggerAnalysis(),
 126   fBackgroundIdentification(0),
 127   fHistBunchCrossing(0),
 128   fHistTriggerPattern(0),
 129   fSkipTriggerClassSelection(0),
 130   fUsingCustomClasses(0),
 131   fSkipV0(0),
 132   fBIFactorA(1),
 133   fBIFactorC(1),
 134   fBIFactorAC(1),
 135   fComputeBG(0),
 136   fBGStatOffset(0),
 137   fUseBXNumbers(1),
 138   fUseMuonTriggers(0),
 139   fFillingScheme(""),
 140   fBin0CallBack(""),
 141   fBin0CallBackPointer(0)
 142 {
 143   // constructor
 144
 145   fCollTrigClasses.SetOwner(1);
 146   fBGTrigClasses.SetOwner(1);
 147   fTriggerAnalysis.SetOwner(1);
 148   fHistStatistics[0] = 0;
 149   fHistStatistics[1] = 0;
 150
 151   AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kWarning);
 152 }
 153
 154 AliPhysicsSelection::~AliPhysicsSelection()
 155 {
 156   // destructor
 157
 158   fCollTrigClasses.Delete();
 159   fBGTrigClasses.Delete();
 160   fTriggerAnalysis.Delete();
 161
 162   if (fHistStatistics[0])
 163   {
 164     delete fHistStatistics[0];
 165     fHistStatistics[0] = 0;
 166   }
 167   if (fHistStatistics[1])
 168   {
 169     delete fHistStatistics[1];
 170     fHistStatistics[1] = 0;
 171   }
 172
 173   if (fHistBunchCrossing)
 174   {
 175     delete fHistBunchCrossing;
 176     fHistBunchCrossing = 0;
 177   }
 178   if (fHistTriggerPattern)
 179   {
 180     delete fHistTriggerPattern;
 181     fHistTriggerPattern = 0;
 182   }
 183   if (fBackgroundIdentification)
 184   {
 185     delete fBackgroundIdentification;
 186     fBackgroundIdentification = 0;
 187   }
 188 }
 189
 190 UInt_t AliPhysicsSelection::CheckTriggerClass(const AliESDEvent* aEsd, const char* trigger, Int_t& triggerLogic) const
 191 {
 192   // checks if the given trigger class(es) are found for the current event
 193   // format of trigger: +TRIGGER1 -TRIGGER2 [#XXX] [&YY] [*ZZ]
 194   //   requires TRIGGER1 and rejects TRIGGER2
 195   //   in bunch crossing XXX
 196   //   if successful, YY is returned (for association between entry in fCollTrigClasses and AliVEvent::EOfflineTriggerTypes)
 197   //   triggerLogic is filled with ZZ, defaults to kCINT1
 198
 199   Bool_t foundBCRequirement = kFALSE;
 200   Bool_t foundCorrectBC = kFALSE;
 201
 202   UInt_t returnCode = AliVEvent::kUserDefined;
 203   triggerLogic = kCINT1;
 204
 205   TString str(trigger);
 206   TObjArray* tokens = str.Tokenize(" ");
 207
 208   for (Int_t i=0; i < tokens->GetEntries(); i++)
 209   {
 210     TString str2(((TObjString*) tokens->At(i))->String());
 211
 212     if (str2[0] == '+' || str2[0] == '-')
 213     {
 214       Bool_t flag = (str2[0] == '+');
 215
 216       str2.Remove(0, 1);
 217
 218       if (flag && !aEsd->IsTriggerClassFired(str2))
 219       {
 220         AliDebug(AliLog::kDebug, Form("Rejecting event because trigger class %s is not present", str2.Data()));
 221         delete tokens;
 222         return kFALSE;
 223       }
 224       if (!flag && aEsd->IsTriggerClassFired(str2))
 225       {
 226         AliDebug(AliLog::kDebug, Form("Rejecting event because trigger class %s is present", str2.Data()));
 227         delete tokens;
 228         return kFALSE;
 229       }
 230     }
 231     else if (str2[0] == '#')
 232     {
 233       foundBCRequirement = kTRUE;
 234
 235       str2.Remove(0, 1);
 236
 237       Int_t bcNumber = str2.Atoi();
 238       AliDebug(AliLog::kDebug, Form("Checking for bunch crossing number %d", bcNumber));
 239
 240       if (aEsd->GetBunchCrossNumber() == bcNumber)
 241       {
 242         foundCorrectBC = kTRUE;
 243         AliDebug(AliLog::kDebug, Form("Found correct bunch crossing %d", bcNumber));
 244       }
 245     }
 246     else if (str2[0] == '&' && !fUsingCustomClasses)
 247     {
 248       str2.Remove(0, 1);
 249
 250       returnCode = str2.Atoll();
 251     }
 252     else if (str2[0] == '*')
 253     {
 254       str2.Remove(0, 1);
 255
 256       triggerLogic = str2.Atoi();
 257     }
 258     else
 259       AliFatal(Form("Invalid trigger syntax: %s", trigger));
 260   }
 261
 262   delete tokens;
 263
 264   if (foundBCRequirement && !foundCorrectBC)
 265     return kFALSE;
 266
 267   return returnCode;
 268 }
 269
 270 //______________________________________________________________________________
 271 TObject *AliPhysicsSelection::GetStatistics(Option_t *option) const
 272 {
 273 // Get the statistics histograms ("ALL" and "BIN0")
 274    TString opt(option);
 275    opt.ToUpper();
 276    Int_t ihist = 0;
 277    if (opt == "ALL") ihist = kStatIdxAll;
 278    if (opt == "BIN0") ihist = kStatIdxBin0;
 279    return fHistStatistics[ihist];
 280 }
 281
 282 //______________________________________________________________________________
 283 UInt_t AliPhysicsSelection::IsCollisionCandidate(const AliESDEvent* aEsd)
 284 {
 285   // checks if the given event is a collision candidate
 286   //
 287   // returns a bit word describing the fired offline triggers (see AliVEvent::EOfflineTriggerTypes)
 288
 289   if (fCurrentRun != aEsd->GetRunNumber()) {
 290     if (!Initialize(aEsd))
 291       AliFatal(Form("Could not initialize for run %d", aEsd->GetRunNumber()));
 292     if(fComputeBG) SetBIFactors(aEsd); // is this safe here?
 293   }
 294   const AliESDHeader* esdHeader = aEsd->GetHeader();
 295   if (!esdHeader)
 296   {
 297     AliError("ESD Header could not be retrieved");
 298     return kFALSE;
 299   }
 300
 301   // check event type; should be PHYSICS = 7 for data and 0 for MC
 302   if (!fMC)
 303   {
 304     if (esdHeader->GetEventType() != 7)
 305       return kFALSE;
 306   }
 307   else
 308   {
 309     if (esdHeader->GetEventType() != 0)
 310       AliFatal(Form("Invalid event type for MC: %d", esdHeader->GetEventType()));
 311   }
 312
 313   UInt_t accept = 0;
 314
 315   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
 316   for (Int_t i=0; i < count; i++)
 317   {
 318     const char* triggerClass = 0;
 319     if (i < fCollTrigClasses.GetEntries())
 320       triggerClass = ((TObjString*) fCollTrigClasses.At(i))->String();
 321     else
 322       triggerClass = ((TObjString*) fBGTrigClasses.At(i - fCollTrigClasses.GetEntries()))->String();
 323
 324     AliDebug(AliLog::kDebug, Form("Processing trigger class %s", triggerClass));
 325
 326     AliTriggerAnalysis* triggerAnalysis = static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i));
 327
 328     triggerAnalysis->FillTriggerClasses(aEsd);
 329
 330     Int_t triggerLogic = 0;
 331     UInt_t singleTriggerResult = CheckTriggerClass(aEsd, triggerClass, triggerLogic);
 332
 333     if (singleTriggerResult)
 334     {
 335       triggerAnalysis->FillHistograms(aEsd);
 336
 337       Bool_t isBin0 = kFALSE;
 338       if (fBin0CallBack != "") {
 339         AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 340         if (!mgr) {
 341           AliError("Cannot get the analysis manager");
 342         }
 343         else {
 344           isBin0 = ((AliAnalysisTaskSE*)mgr->GetTask(fBin0CallBack.Data()))->IsEventInBinZero();
 345         }
 346       } else if (fBin0CallBackPointer) {
 347           isBin0 = (*fBin0CallBackPointer)(aEsd);
 348       }
 349
 350       // hardware trigger
 351       Int_t fastORHW = triggerAnalysis->SPDFiredChips(aEsd, 1); // SPD number of chips from trigger bits (!)
 352       Int_t fastORHWL1 = triggerAnalysis->SPDFiredChips(aEsd, 1, kFALSE, 2); // SPD number of chips from trigger bits in second layer (!)
 353       Bool_t v0AHW     = fSkipV0 ? 0 :(triggerAnalysis->V0Trigger(aEsd, AliTriggerAnalysis::kASide, kTRUE) == AliTriggerAnalysis::kV0BB);// should replay hw trigger
 354       Bool_t v0CHW     = fSkipV0 ? 0 :(triggerAnalysis->V0Trigger(aEsd, AliTriggerAnalysis::kCSide, kTRUE) == AliTriggerAnalysis::kV0BB);// should replay hw trigger
 355
 356       // offline trigger
 357       Int_t fastOROffline = triggerAnalysis->SPDFiredChips(aEsd, 0); // SPD number of chips from clusters (!)
 358       Int_t fastOROfflineL1 = triggerAnalysis->SPDFiredChips(aEsd, 0, kFALSE, 2); // SPD number of chips from clusters in second layer (!)
 359       Bool_t v0A       = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0A);
 360       Bool_t v0C       = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0C);
 361       Bool_t v0ABG = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0ABG);
 362       Bool_t v0CBG = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0CBG);
 363       Bool_t v0BG = v0ABG || v0CBG;
 364
 365       // fmd
 366       Bool_t fmdA =  triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kFMDA);
 367       Bool_t fmdC =  triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kFMDC);
 368       Bool_t fmd  = fmdA || fmdC;
 369
 370       // SSD
 371       //Int_t ssdClusters = triggerAnalysis->SSDClusters(aEsd);
 372
 373       // ZDC
 374       Bool_t zdcA = triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kZDCA);
 375       Bool_t zdcC = triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kZDCC);
 376
 377
 378       // Some "macros"
 379       Bool_t mb1 = (fastOROffline > 0 || v0A || v0C) && (!v0BG);
 380       Bool_t mb1prime = (fastOROffline > 1 || (fastOROffline > 0 && (v0A || v0C)) || (v0A && v0C) ) && (!v0BG);
 381
 382       // Background rejection
 383       Bool_t bgID = kFALSE;
 384       if (fBackgroundIdentification)
 385         bgID = ! fBackgroundIdentification->IsSelected(const_cast<AliESDEvent*> (aEsd));
 386
 387       /*Int_t ntrig = fastOROffline; // any 2 hits
 388       if(v0A)              ntrig += 1;
 389       if(v0C)              ntrig += 1; //v0C alone is enough
 390       if(fmd)              ntrig += 1;
 391       if(ssdClusters>1)    ntrig += 1;*/
 392
 393       // replay hardware trigger (should only remove events for MC)
 394
 395       Bool_t hwTrig = kFALSE;
 396       switch (triggerLogic)
 397       {
 398         case kCINT1:  hwTrig = fastORHW > 0 || v0AHW || v0CHW; break;
 399         case kCMBS2A: hwTrig = fastORHWL1 > 1 && v0AHW; break;
 400         case kCMBS2C: hwTrig = fastORHWL1 > 1 && v0CHW; break;
 401         case kCMBAC:  hwTrig = v0AHW && v0CHW; break;
 402         default: AliFatal(Form("Undefined trigger logic %d", triggerLogic)); break;
 403       }
 404
 405       // Fill trigger pattern histo
 406       Int_t tpatt = 0;
 407       if (fastORHW>0) tpatt+=1;
 408       if (v0AHW)      tpatt+=2;
 409       if (v0CHW)      tpatt+=4;
 410       fHistTriggerPattern->Fill( tpatt );
 411
 412       // fill statistics and return decision
 413       const Int_t nHistStat = 2;
 414       for(Int_t iHistStat = 0; iHistStat < nHistStat; iHistStat++){
 415         if (iHistStat == kStatIdxBin0 && !isBin0) continue; // skip the filling of bin0 stats if the event is not in the bin0
 416
 417         fHistStatistics[iHistStat]->Fill(kStatTriggerClass, i);
 418
 419         // We fill the rest only if hw trigger is ok
 420         if (!hwTrig)
 421           {
 422             AliDebug(AliLog::kDebug, "Rejecting event because hardware trigger is not fired");
 423             continue;
 424           } else {
 425           fHistStatistics[iHistStat]->Fill(kStatHWTrig, i);
 426         }
 427
 428
 429         // v0 BG stats
 430         if (v0ABG)
 431           fHistStatistics[iHistStat]->Fill(kStatV0ABG, i);
 432         if (v0CBG)
 433           fHistStatistics[iHistStat]->Fill(kStatV0CBG, i);
 434
 435         // We fill the rest only if mb1 && ! v0BG
 436         if (mb1)
 437           fHistStatistics[iHistStat]->Fill(kStatMB1, i);
 438         else continue;
 439
 440         if (mb1prime)
 441           fHistStatistics[iHistStat]->Fill(kStatMB1Prime, i);
 442
 443         if (fmd)
 444           fHistStatistics[iHistStat]->Fill(kStatFMD, i);
 445
 446         //if(ntrig >= 2 && !v0BG)
 447         //  fHistStatistics[iHistStat]->Fill(kStatAny2Hits, i);
 448
 449         if (fastOROffline > 0)
 450           fHistStatistics[iHistStat]->Fill(kStatFO1, i);
 451         if (fastOROffline > 1)
 452           fHistStatistics[iHistStat]->Fill(kStatFO2, i);
 453         if (fastOROfflineL1 > 1)
 454           fHistStatistics[iHistStat]->Fill(kStatFO2L1, i);
 455
 456         if (v0A)
 457           fHistStatistics[iHistStat]->Fill(kStatV0A, i);
 458         if (v0C)
 459           fHistStatistics[iHistStat]->Fill(kStatV0C, i);
 460
 461         if (zdcA)
 462           fHistStatistics[iHistStat]->Fill(kStatZDCA, i);
 463         if (zdcC)
 464           fHistStatistics[iHistStat]->Fill(kStatZDCC, i);
 465         if (zdcA && zdcC)
 466           fHistStatistics[iHistStat]->Fill(kStatZDCAC, i);
 467
 468         //       if (fastOROffline > 1 && !v0BG)
 469         //         fHistStatistics[iHistStat]->Fill(kStatFO2NoBG, i);
 470
 471         //if (fastOROffline > 0 && (v0A || v0C) && !v0BG)
 472         //  fHistStatistics[iHistStat]->Fill(kStatFO1AndV0, i);
 473
 474         if (v0A && v0C && !v0BG && !bgID)
 475           fHistStatistics[iHistStat]->Fill(kStatV0, i);
 476
 477         Bool_t offlineAccepted = kFALSE;
 478
 479         switch (triggerLogic)
 480         {
 481           case kCINT1:  offlineAccepted = mb1; break;
 482           case kCMBS2A: offlineAccepted = fastOROfflineL1 > 1 && v0A; break;
 483           case kCMBS2C: offlineAccepted = fastOROfflineL1 > 1 && v0C; break;
 484           case kCMBAC:  offlineAccepted = v0A && v0C; break;
 485           default: AliFatal(Form("Undefined trigger logic %d", triggerLogic)); break;
 486         }
 487
 488         if ( offlineAccepted )
 489           {
 490             if (!v0BG || fSkipV0)
 491               {
 492                 if (!v0BG) fHistStatistics[iHistStat]->Fill(kStatOffline, i);
 493
 494                 if (fBackgroundIdentification && bgID)
 495                   {
 496                     AliDebug(AliLog::kDebug, "Rejecting event because of background identification");
 497                     fHistStatistics[iHistStat]->Fill(kStatBG, i);
 498                   }
 499                 else
 500                   {
 501                     AliDebug(AliLog::kDebug, Form("Accepted event for histograms with trigger logic %d", triggerLogic));
 502
 503                     fHistStatistics[iHistStat]->Fill(kStatAccepted, i);
 504                     if(iHistStat == kStatIdxAll) fHistBunchCrossing->Fill(aEsd->GetBunchCrossNumber(), i); // Fill only for all (avoid double counting)
 505                     if((i < fCollTrigClasses.GetEntries() || fSkipTriggerClassSelection) && (iHistStat==kStatIdxAll))
 506                       accept |= singleTriggerResult; // only set for "all" (should not really matter)
 507                   }
 508               }
 509             else
 510               AliDebug(AliLog::kDebug, "Rejecting event because of V0 BG flag");
 511           }
 512         else
 513           AliDebug(AliLog::kDebug, Form("Rejecting event because trigger logic %d is not fulfilled", triggerLogic));
 514       }
 515     }
 516   }
 517
 518   if (accept)
 519     AliDebug(AliLog::kDebug, Form("Accepted event as collision candidate with bit mask %d", accept));
 520
 521   return accept;
 522 }
 523
 524 Int_t AliPhysicsSelection::GetTriggerScheme(UInt_t runNumber) const
 525 {
 526   // returns the current trigger scheme (classes that are accepted/rejected)
 527
 528   if (fMC)
 529     return 0;
 530
 531   // TODO dependent on run number
 532
 533   switch (runNumber)
 534   {
 535     // CSMBB triggers
 536     case 104044:
 537     case 105054:
 538     case 105057:
 539       return 2;
 540   }
 541
 542   // defaults
 543   return 1;
 544 }
 545
 546 const char * AliPhysicsSelection::GetFillingScheme(UInt_t runNumber)  {
 547
 548   if(fMC) return "MC";
 549
 550   if      (runNumber >= 104065 && runNumber <= 104160) {
 551     return "4x4a";
 552   }
 553   else if (runNumber >= 104315 && runNumber <= 104321) {
 554     return "4x4a*";
 555   }
 556   else if (runNumber >= 104792 && runNumber <= 104803) {
 557     return "4x4b";
 558   }
 559   else if (runNumber >= 104824 && runNumber <= 104892) {
 560     return "4x4c";
 561   }
 562   else if (runNumber == 105143 || runNumber == 105160) {
 563     return "16x16a";
 564   }
 565   else if (runNumber >= 105256 && runNumber <= 105268) {
 566     return "4x4c";
 567   }
 568   else if (runNumber >= 114786 && runNumber <= 116684) {
 569     return "Single_2b_1_1_1";
 570   }
 571   else if (runNumber >= 117048 && runNumber <= 117120) {
 572     return "Single_3b_2_2_2";
 573   }
 574   else if (runNumber >= 117220 && runNumber <= 119163) {
 575     return "Single_2b_1_1_1";
 576   }
 577   else if (runNumber >= 119837 && runNumber <= 119862) {
 578     return "Single_4b_2_2_2";
 579   }
 580   else if (runNumber >= 119902 && runNumber <= 120691) {
 581     return "Single_6b_3_3_3";
 582   }
 583   else if (runNumber >= 120741 && runNumber <= 122375) {
 584     return "Single_13b_8_8_8";
 585   }
 586   else if (runNumber >= 130148 && runNumber <= 130375) {
 587     return "125n_48b_36_16_36";
 588   }
 589   else if (runNumber >= 130601 && runNumber <= 130640) {
 590     return "1000ns_50b_35_14_35";
 591   }
 592   else {
 593     AliError(Form("Unknown filling scheme (run %d)", runNumber));
 594   }
 595
 596   return "Unknown";
 597 }
 598
 599 const char * AliPhysicsSelection::GetBXIDs(UInt_t runNumber, const char * trigger)  {
 600
 601   if (!fUseBXNumbers || fMC) return "";
 602
 603   if      (runNumber >= 104065 && runNumber <= 104160) {
 604     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #2128 #3019";
 605     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #346 #3465";
 606     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1234 #1680";
 607     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
 608     //    else AliError(Form("Unknown trigger: %s", trigger));
 609   }
 610   else if (runNumber >= 104315 && runNumber <= 104321) {
 611     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #2000 #2891";
 612     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #218 #3337";
 613     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1106 #1552";
 614     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
 615     //    else AliError(Form("Unknown trigger: %s", trigger));
 616   }
 617   else if (runNumber >= 104792 && runNumber <= 104803) {
 618     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #2228 #3119";
 619     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2554 #446";
 620     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1334 #769";
 621     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
 622     //    else AliError(Form("Unknown trigger: %s", trigger));
 623   }
 624   else if (runNumber >= 104824 && runNumber <= 104892) {
 625     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #3119 #769";
 626     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2554 #446";
 627     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1334 #2228";
 628     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
 629     //    else AliError(Form("Unknown trigger: %s", trigger));
 630   }
 631   else if (runNumber == 105143 || runNumber == 105160) {
 632     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #1337 #1418 #2228 #2309 #3119 #3200 #446 #527";
 633     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #1580  #1742  #1904  #2066  #2630  #2792  #2954  #3362";
 634     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "  #845  #1007  #1169   #1577 #3359 #3521 #119  #281 ";
 635     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
 636     //    else AliError(Form("Unknown trigger: %s", trigger));
 637   }
 638   else if (runNumber >= 105256 && runNumber <= 105268) {
 639     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #3019 #669";
 640     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2454 #346";
 641     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1234 #2128";
 642     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
 643     //    else AliError(Form("Unknown trigger: %s", trigger));
 644   } else if (runNumber >= 114786 && runNumber <= 116684) { // 7 TeV 2010, assume always the same filling scheme
 645     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #346";
 646     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2131";
 647     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #3019";
 648     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #1238";
 649     //    else AliError(Form("Unknown trigger: %s", trigger));
 650   }
 651   else if (runNumber >= 117048 && runNumber <= 117120) {
 652     //    return "Single_3b_2_2_2";
 653    if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346  #1240 ";
 654    else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #2131 ";
 655    else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019 ";
 656    else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1238";
 657    //   else AliError(Form("Unknown trigger: %s", trigger));
 658
 659   }
 660   else if ((runNumber >= 117220 && runNumber <= 118555) || (runNumber >= 118784 && runNumber <= 119163))
 661   {
 662     //    return "Single_2b_1_1_1";
 663     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346 ";
 664     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #2131 ";
 665     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019 ";
 666     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1238 ";
 667     //    else AliError(Form("Unknown trigger: %s", trigger));
 668   }
 669   else if (runNumber >= 118556 && runNumber <= 118783) {
 670     //    return "Single_2b_1_1_1";
 671     // same as previous but was misaligned by 1 BX in fill 1069
 672     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #345 ";
 673     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #2130 ";
 674     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3018 ";
 675     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1238 ";
 676     //    else AliError(Form("Unknown trigger: %s", trigger));
 677   }
 678   else if (runNumber >= 119837 && runNumber <= 119862) {
 679     //    return "Single_4b_2_2_2";
 680     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #669  #3019 ";
 681     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #346  #2454 ";
 682     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #1234  #2128 ";
 683     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1681 #3463";
 684     //    else AliError(Form("Unknown trigger: %s", trigger));
 685
 686   }
 687   else if (runNumber >= 119902 && runNumber <= 120691) {
 688     //    return "Single_6b_3_3_3";
 689     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346  #546  #746 ";
 690     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #2131  #2331  #2531 ";
 691     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019  #3219  #3419 ";
 692     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1296 #1670";
 693     //    else AliError(Form("Unknown trigger: %s", trigger));
 694   }
 695   else if (runNumber >= 120741 && runNumber <= 122375) {
 696     //    return "Single_13b_8_8_8";
 697     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346  #446  #546  #646  #1240  #1340  #1440  #1540 ";
 698     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #946  #2131  #2231  #2331  #2431 ";
 699     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019  #3119  #3219  #3319  #3519 ";
 700     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1835 #2726";
 701     //    else AliError(Form("Unknown trigger: %s", trigger));
 702
 703   }
 704   else if (runNumber >= 130148 && runNumber <= 130375) {
 705     TString triggerString = trigger;
 706     static TString returnString = " ";
 707     returnString = "";
 708     if (triggerString.Contains("B")) returnString += "   #346  #396  #446  #496  #546  #596  #646  #696  #1240  #1290  #1340  #1390  #1440  #1490  #1540  #1590 ";
 709     if (triggerString.Contains("A")) returnString += "   #755  #805  #855  #905  #955  #1005  #1799  #1849  #1899  #2131  #2181  #2231  #2281  #2331  #2381  #2431  #2481  #2531  #2581  #2631  #2846  #3016  #3066  #3116  #3166  #3216  #3266  #3316  #3366  #3425  #3475  #3525 ";
 710     if (triggerString.Contains("C")) returnString += "   #3019  #3069  #3119  #3169  #3219  #3269  #3319  #3369  #14  #64  #114  #746  #796  #846  #908  #958  #1008  #1640  #1690  #1740  #2055  #2125  #2175  #2225  #2275  #2325  #2375  #2425  #2475  #2534  #2584  #2634 ";
 711     // Printf("0x%x",returnString.Data());
 712     // Printf("%s",returnString.Data());
 713     return returnString.Data();
 714   }
 715   else if (runNumber >= 130601 && runNumber <= 130640) {
 716     TString triggerString = trigger;
 717     static TString returnString = " ";
 718     returnString = "";
 719     if (triggerString.Contains("B")) returnString += "  #346  #386  #426  #466  #506  #546  #586  #1240  #1280  #1320  #1360  #1400  #1440  #1480 ";
 720     if (triggerString.Contains("A")) returnString += "  #626  #666  #706  #746  #786  #826  #866  #1520  #1560  #1600  #1640  #1680  #1720  #1760  #2076  #2131  #2171  #2211  #2251  #2291  #2331  #2371  #2414  #2454  #2494  #2534  #2574  #2614  #2654  #2694  #2734  #2774  #2814 "; //#2854  #2894  #2934 not present in this run
 721     if (triggerString.Contains("C")) returnString += "  #3019  #3059  #3099  #3139  #3179  #3219  #3259  #3299  #3339  #3379  #3419  #3459  #3499  #3539  #115  #629  #669  #709  #749  #789  #829  #869  #909  #949  #989  #1029  #1069  #1109  #1149  #1523  #1563  #1603  #1643 "; //#1683  #1723  #1763 not present in this run
 722     return returnString.Data();
 723   }
 724
 725   else {
 726     AliWarning(Form("Unknown run %d, using all BXs!",runNumber));
 727   }
 728
 729   return "";
 730 }
 731
 732 Bool_t AliPhysicsSelection::Initialize(const AliESDEvent* aEsd)
 733 {
 734   // initializes the object for the given ESD
 735
 736   AliInfo(Form("Initializing for beam type: %s", aEsd->GetESDRun()->GetBeamType()));
 737   Bool_t pp = kTRUE;
 738   if (strcmp(aEsd->GetESDRun()->GetBeamType(), "Pb-Pb") == 0)
 739     pp = kFALSE;
 740
 741   return Initialize(aEsd->GetRunNumber(), pp);
 742 }
 743
 744 Bool_t AliPhysicsSelection::Initialize(Int_t runNumber, Bool_t pp)
 745 {
 746   // initializes the object for the given run
 747   // TODO having the run number here and parameters hardcoded is clearly temporary, a way needs to be found to have a CDB-like configuration also for analysis
 748
 749   Bool_t oldStatus = TH1::AddDirectoryStatus();
 750   TH1::AddDirectory(kFALSE);
 751
 752   if(!fBin0CallBack)
 753     AliError("Bin0 Callback not set: will not fill the statistics for the bin 0");
 754
 755   if (fMC) {
 756     // override BX and bg options in case of MC
 757     fComputeBG    = kFALSE;
 758     fUseBXNumbers = kFALSE;
 759   }
 760
 761   Int_t triggerScheme = GetTriggerScheme(runNumber);
 762   if (!fUsingCustomClasses && fCurrentRun != -1 && triggerScheme != GetTriggerScheme(fCurrentRun))
 763     AliFatal("Processing several runs with different trigger schemes is not supported");
 764
 765   if(fComputeBG && fCurrentRun != -1 && fCurrentRun != runNumber)
 766     AliFatal("Cannot process several runs because BG computation is requested");
 767
 768   if(fComputeBG && !fUseBXNumbers)
 769     AliFatal("Cannot compute BG if BX numbers are not used");
 770
 771   if(fUseBXNumbers && fFillingScheme != "" && fFillingScheme != GetFillingScheme(runNumber))
 772     AliFatal("Cannot process runs with different filling scheme if usage of BX numbers is requested");
 773
 774   fFillingScheme      = GetFillingScheme(runNumber);
 775
 776   AliInfo(Form("Initializing for run %d", runNumber));
 777
 778   // initialize first time?
 779   if (fCurrentRun == -1)
 780   {
 781     if (fUsingCustomClasses) {
 782       AliInfo("Using user-provided trigger classes");
 783     } else {
 784       if (pp)
 785       {
 786         switch (triggerScheme)
 787         {
 788         case 0:
 789           // MC Proton-Proton
 790
 791           fCollTrigClasses.Add(new TObjString(Form("&%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCINT1)));
 792           break;
 793
 794         case 1:
 795           // Proton-Proton
 796
 797           // trigger classes used before August 2010
 798           fCollTrigClasses.Add(new TObjString(Form("%s%s &%u","+CINT1B-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1B-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
 799           fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CINT1A-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1A-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
 800           fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CINT1C-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1C-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
 801           fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CINT1-E-NOPF-ALL",      GetBXIDs(runNumber,"CINT1-E-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
 802
 803           // Muon trigger have the same BXIDs of the corresponding CINT triggers
 804           fCollTrigClasses.Add(new TObjString(Form("%s%s &%u","+CMUS1B-ABCE-NOPF-MUON",  GetBXIDs(runNumber,"CINT1B-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
 805           fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CMUS1A-ABCE-NOPF-MUON",  GetBXIDs(runNumber,"CINT1A-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
 806           fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CMUS1C-ABCE-NOPF-MUON",  GetBXIDs(runNumber,"CINT1C-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
 807           fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CMUS1-E-NOPF-MUON"    ,  GetBXIDs(runNumber,"CINT1-E-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
 808
 809           // triggers classes used from August 2010
 810           // MB
 811           fCollTrigClasses.Add(new TObjString(Form("%s%s &%u", "+CINT1-B-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "B"),  (UInt_t) AliVEvent::kMB)));
 812           fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CINT1-AC-NOPF-ALLNOTRD", GetBXIDs(runNumber, "AC"), (UInt_t) AliVEvent::kMB)));
 813           fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CINT1-E-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "E"),  (UInt_t) AliVEvent::kMB)));
 814
 815           // MUON
 816           fCollTrigClasses.Add(new TObjString(Form("%s%s &%u", "+CMUS1-B-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "B"),  (UInt_t) AliVEvent::kMUON)));
 817           fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CMUS1-AC-NOPF-ALLNOTRD", GetBXIDs(runNumber, "AC"), (UInt_t) AliVEvent::kMUON)));
 818           fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CMUS1-E-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "E"),  (UInt_t) AliVEvent::kMUON)));
 819
 820           // High Multiplicity
 821           fCollTrigClasses.Add(new TObjString(Form("%s%s &%u", "+CSH1-B-NOPF-ALLNOTRD"  , GetBXIDs(runNumber, "B"),  (UInt_t) AliVEvent::kHighMult)));
 822           fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CSH1-AC-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "AC"), (UInt_t) AliVEvent::kHighMult)));
 823           fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CSH1-E-NOPF-ALLNOTRD"  , GetBXIDs(runNumber, "E"),  (UInt_t) AliVEvent::kHighMult)));
 824
 825           // WARNING: IF YOU ADD MORE TRIGGER CLASSES, PLEASE CHECK THAT THE REGULAR EXPRESSION IN GetStatRow IS STILL VALID
 826
 827           break;
 828
 829         case 2:
 830           // Proton-Proton
 831
 832           fCollTrigClasses.Add(new TObjString(Form("+CSMBB-ABCE-NOPF-ALL &%u", (UInt_t) AliVEvent::kMB)));
 833           fBGTrigClasses.Add(new TObjString(Form("+CSMBA-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL &%u", (UInt_t) AliVEvent::kMB)));
 834           fBGTrigClasses.Add(new TObjString(Form("+CSMBC-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL &%u", (UInt_t) AliVEvent::kMB)));
 835           break;
 836
 837         default:
 838           AliFatal(Form("Unsupported trigger scheme %d", triggerScheme));
 839         }
 840       }
 841       else
 842       {
 843         switch (triggerScheme)
 844         {
 845         case 0:
 846           // MC Heavy-Ion
 847
 848           fCollTrigClasses.Add(new TObjString(Form("&%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2A)));
 849           fCollTrigClasses.Add(new TObjString(Form("&%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2C)));
 850           fCollTrigClasses.Add(new TObjString(Form("&%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBAC)));
 851           break;
 852
 853         case 1:
 854           // Data Heavy-ion
 855
 856           fCollTrigClasses.Add(new TObjString(Form("+CMBS2A-B-NOPF-ALL &%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2A)));
 857           fBGTrigClasses.Add  (new TObjString(Form("+CMBS2A-A-NOPF-ALL &%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2A)));
 858           fBGTrigClasses.Add  (new TObjString(Form("+CMBS2A-C-NOPF-ALL &%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2A)));
 859           fBGTrigClasses.Add  (new TObjString(Form("+CMBS2A-E-NOPF-ALL &%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2A)));
 860
 861           fCollTrigClasses.Add(new TObjString(Form("+CMBS2C-B-NOPF-ALL &%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2C)));
 862           fBGTrigClasses.Add  (new TObjString(Form("+CMBS2C-A-NOPF-ALL &%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2C)));
 863           fBGTrigClasses.Add  (new TObjString(Form("+CMBS2C-C-NOPF-ALL &%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2C)));
 864           fBGTrigClasses.Add  (new TObjString(Form("+CMBS2C-E-NOPF-ALL &%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2C)));
 865
 866           fCollTrigClasses.Add(new TObjString(Form("+CMBAC-B-NOPF-ALL &%u *%d",  (UInt_t) AliVEvent::kMB, (Int_t) kCMBAC)));
 867           fBGTrigClasses.Add  (new TObjString(Form("+CMBAC-A-NOPF-ALL &%u *%d",  (UInt_t) AliVEvent::kMB, (Int_t) kCMBAC)));
 868           fBGTrigClasses.Add  (new TObjString(Form("+CMBAC-C-NOPF-ALL &%u *%d",  (UInt_t) AliVEvent::kMB, (Int_t) kCMBAC)));
 869           fBGTrigClasses.Add  (new TObjString(Form("+CMBAC-E-NOPF-ALL &%u *%d",  (UInt_t) AliVEvent::kMB, (Int_t) kCMBAC)));
 870
 871           break;
 872
 873         default:
 874           AliFatal(Form("Unsupported trigger scheme %d", triggerScheme));
 875         }
 876       }
 877     }
 878
 879     Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
 880
 881     for (Int_t i=0; i<count; i++)
 882     {
 883       AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis;
 884       triggerAnalysis->SetAnalyzeMC(fMC);
 885       triggerAnalysis->EnableHistograms();
 886       triggerAnalysis->SetSPDGFOThreshhold(1);
 887       triggerAnalysis->SetDoFMD(kFALSE);
 888       fTriggerAnalysis.Add(triggerAnalysis);
 889     }
 890
 891     // TODO: shall I really delete this?
 892     if (fHistStatistics[0])
 893       delete fHistStatistics[0];
 894     if (fHistStatistics[1])
 895       delete fHistStatistics[1];
 896
 897     fHistStatistics[kStatIdxBin0] = BookHistStatistics("_Bin0");
 898     fHistStatistics[kStatIdxAll]  = BookHistStatistics("");
 899
 900     if (fHistBunchCrossing)
 901       delete fHistBunchCrossing;
 902
 903     fHistBunchCrossing = new TH2F("fHistBunchCrossing", "fHistBunchCrossing;bunch crossing number;", 4000, -0.5, 3999.5,  count, -0.5, -0.5 + count);
 904
 905     if (fHistTriggerPattern)
 906       delete fHistTriggerPattern;
 907
 908     const int ntrig=3;
 909     Int_t n = 1;
 910     const Int_t nbinTrig = TMath::Nint(TMath::Power(2,ntrig));
 911
 912     fHistTriggerPattern = new TH1F("fHistTriggerPattern", "Trigger pattern: FO + 2*v0A + 4*v0C",
 913                                    nbinTrig, -0.5, nbinTrig-0.5);
 914     fHistTriggerPattern->GetXaxis()->SetBinLabel(1,"NO TRIG");
 915     fHistTriggerPattern->GetXaxis()->SetBinLabel(2,"FO");
 916     fHistTriggerPattern->GetXaxis()->SetBinLabel(3,"v0A");
 917     fHistTriggerPattern->GetXaxis()->SetBinLabel(4,"FO & v0A");
 918     fHistTriggerPattern->GetXaxis()->SetBinLabel(5,"v0C");
 919     fHistTriggerPattern->GetXaxis()->SetBinLabel(6,"FO & v0C");
 920     fHistTriggerPattern->GetXaxis()->SetBinLabel(7,"v0A & v0C");
 921     fHistTriggerPattern->GetXaxis()->SetBinLabel(8,"FO & v0A & v0C");
 922
 923
 924     n = 1;
 925     for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
 926     {
 927       fHistBunchCrossing->GetYaxis()->SetBinLabel(n, ((TObjString*) fCollTrigClasses.At(i))->String());
 928       n++;
 929     }
 930     for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
 931     {
 932       fHistBunchCrossing->GetYaxis()->SetBinLabel(n, ((TObjString*) fBGTrigClasses.At(i))->String());
 933       n++;
 934     }
 935
 936
 937
 938   }
 939
 940   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
 941   for (Int_t i=0; i<count; i++)
 942   {
 943     AliTriggerAnalysis* triggerAnalysis = static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i));
 944
 945     switch (runNumber)
 946     {
 947       case 104315:
 948       case 104316:
 949       case 104320:
 950       case 104321:
 951       case 104439:
 952         triggerAnalysis->SetV0TimeOffset(7.5);
 953         break;
 954       default:
 955         triggerAnalysis->SetV0TimeOffset(0);
 956     }
 957   }
 958
 959   fCurrentRun = runNumber;
 960
 961   TH1::AddDirectory(oldStatus);
 962
 963   return kTRUE;
 964 }
 965
 966 TH2F * AliPhysicsSelection::BookHistStatistics(const char * tag) {
 967   // add 6 rows to count for the estimate of good, accidentals and
 968   // BG and the ratio of BG and accidentals to total +ratio goot to
 969   // first col + 2 for error on good.
 970   // TODO: Remember the the indexes of rows for the BG selection. Add new member fBGRows[] and use kStat as indexes
 971
 972   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
 973 #ifdef VERBOSE_STAT
 974   Int_t extrarows = fComputeBG != 0 ? 8 : 0;
 975 #else
 976   Int_t extrarows = fComputeBG != 0 ? 3 : 0;
 977 #endif
 978   TH2F * h = new TH2F(Form("fHistStatistics%s",tag), Form("fHistStatistics - %s ;;",tag), kStatAccepted, 0.5, kStatAccepted+0.5, count+extrarows, -0.5, -0.5 + count+extrarows);
 979
 980   h->GetXaxis()->SetBinLabel(kStatTriggerClass,  "Trigger class");
 981   h->GetXaxis()->SetBinLabel(kStatHWTrig,        "Hardware trigger");
 982   h->GetXaxis()->SetBinLabel(kStatFO1,           "FO >= 1");
 983   h->GetXaxis()->SetBinLabel(kStatFO2,           "FO >= 2");
 984   h->GetXaxis()->SetBinLabel(kStatFO2L1,         "FO (L1) >= 2");
 985   h->GetXaxis()->SetBinLabel(kStatV0A,           "V0A");
 986   h->GetXaxis()->SetBinLabel(kStatV0C,           "V0C");
 987   h->GetXaxis()->SetBinLabel(kStatFMD,           "FMD");
 988   h->GetXaxis()->SetBinLabel(kStatV0ABG,         "V0A BG");
 989   h->GetXaxis()->SetBinLabel(kStatV0CBG,         "V0C BG");
 990   h->GetXaxis()->SetBinLabel(kStatZDCA,          "ZDCA");
 991   h->GetXaxis()->SetBinLabel(kStatZDCC,          "ZDCC");
 992   h->GetXaxis()->SetBinLabel(kStatZDCAC,         "ZDCA & ZDCC");
 993   h->GetXaxis()->SetBinLabel(kStatMB1,           "(FO >= 1 | V0A | V0C) & !V0 BG");
 994   h->GetXaxis()->SetBinLabel(kStatMB1Prime,      "(FO >= 2 | (FO >= 1 & (V0A | V0C)) | (V0A &v0C) ) & !V0 BG");
 995   //h->GetXaxis()->SetBinLabel(kStatFO1AndV0,    "FO >= 1 & (V0A | V0C) & !V0 BG");
 996   h->GetXaxis()->SetBinLabel(kStatV0,            "V0A & V0C & !V0 BG & !BG ID");
 997   h->GetXaxis()->SetBinLabel(kStatOffline,       "Offline Trigger");
 998   //h->GetXaxis()->SetBinLabel(kStatAny2Hits,    "2 Hits & !V0 BG");
 999   h->GetXaxis()->SetBinLabel(kStatBG,            "Background identification");
1000   h->GetXaxis()->SetBinLabel(kStatAccepted,      "Accepted");
1001
1002   Int_t n = 1;
1003   for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
1004     {
1005       h->GetYaxis()->SetBinLabel(n, ((TObjString*) fCollTrigClasses.At(i))->String());
1006       n++;
1007     }
1008   for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
1009     {
1010       h->GetYaxis()->SetBinLabel(n, ((TObjString*) fBGTrigClasses.At(i))->String());
1011       n++;
1012     }
1013
1014   if(fComputeBG) {
1015     fBGStatOffset = n;
1016     h->GetYaxis()->SetBinLabel(n++, Form("BG (A+C) (Mask [0x%x])", fComputeBG));
1017     h->GetYaxis()->SetBinLabel(n++, "ACC");
1018 #ifdef VERBOSE_STAT
1019     h->GetYaxis()->SetBinLabel(n++, "BG (A+C) %  (rel. to CINT1B)");
1020     h->GetYaxis()->SetBinLabel(n++, "ACC % (rel. to CINT1B)");
1021     h->GetYaxis()->SetBinLabel(n++, "ERR GOOD %");
1022     h->GetYaxis()->SetBinLabel(n++, "GOOD % (rel. to 1st col)");
1023     h->GetYaxis()->SetBinLabel(n++, "ERR GOOD");
1024 #endif
1025     h->GetYaxis()->SetBinLabel(n++, "GOOD");
1026   }
1027
1028   return h;
1029 }
1030
1031 void AliPhysicsSelection::Print(Option_t *option) const
1032 {
1033   // print the configuration
1034   TString msg;
1035   Printf("Configuration initialized for run %d (MC: %d):", fCurrentRun, fMC);
1036   msg += Form("Configuration initialized for run %d (MC: %d):\n", fCurrentRun, fMC);
1037
1038   Printf("Collision trigger classes:");
1039   for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
1040     Printf("%s", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
1041
1042   Printf("Background trigger classes:");
1043   for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
1044     Printf("%s", ((TObjString*) fBGTrigClasses.At(i))->String().Data());
1045
1046   AliTriggerAnalysis* triggerAnalysis = dynamic_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(0));
1047
1048   if (triggerAnalysis)
1049   {
1050     if (triggerAnalysis->GetV0TimeOffset() > 0)
1051       Printf("V0 time offset active: %.2f ns", triggerAnalysis->GetV0TimeOffset());
1052
1053     Printf("\nTotal available events:");
1054
1055     triggerAnalysis->PrintTriggerClasses();
1056   }
1057
1058   if (fHistStatistics[kStatIdxAll])
1059   {
1060     for (Int_t i=0; i<fCollTrigClasses.GetEntries(); i++)
1061     {
1062       Printf("\nSelection statistics for collision trigger %s:", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
1063       msg += Form("\nSelection statistics for collision trigger %s:\n", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
1064
1065       Printf("Total events with correct trigger class: %d", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(1, i+1));
1066       msg += Form("Total events with correct trigger class: %d\n", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(1, i+1));
1067       Printf("Selected collision candidates: %d", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(fHistStatistics[kStatIdxAll]->GetXaxis()->FindBin("Accepted"), i+1));
1068       msg += Form("Selected collision candidates: %d\n", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(fHistStatistics[kStatIdxAll]->GetXaxis()->FindBin("Accepted"), i+1));
1069     }
1070   }
1071
1072   if (fHistBunchCrossing)
1073   {
1074     Printf("\nBunch crossing statistics:");
1075     msg += "\nBunch crossing statistics:\n";
1076
1077     for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
1078     {
1079       TString str;
1080       str.Form("Trigger %s has accepted events in the bunch crossings: ", fHistBunchCrossing->GetYaxis()->GetBinLabel(i));
1081
1082       for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
1083         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
1084           str += Form("%d, ", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
1085
1086       Printf("%s", str.Data());
1087       msg += str;
1088       msg += "\n";
1089     }
1090
1091     for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
1092     {
1093       Int_t countColl = 0;
1094       Int_t countBG = 0;
1095       for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
1096       {
1097         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
1098         {
1099           if (fCollTrigClasses.FindObject(fHistBunchCrossing->GetYaxis()->GetBinLabel(i)))
1100             countColl++;
1101           if (fBGTrigClasses.FindObject(fHistBunchCrossing->GetYaxis()->GetBinLabel(i)))
1102             countBG++;
1103         }
1104       }
1105       if (countColl > 0 && countBG > 0)
1106         Printf("WARNING: Bunch crossing %d has collision and BG trigger classes active. Check BPTX functioning for this run!", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
1107     }
1108   }
1109
1110   if (fUsingCustomClasses)
1111     Printf("WARNING: Using custom trigger classes!");
1112   if (fSkipTriggerClassSelection)
1113     Printf("WARNING: Skipping trigger class selection!");
1114   if (fSkipV0)
1115     Printf("WARNING: Ignoring V0 information in selection");
1116   if(!fBin0CallBack)
1117     Printf("WARNING: Callback not set: will not fill the statistics for the bin 0");
1118   TString opt(option);
1119   opt.ToUpper();
1120   if (opt == "STAT") {
1121      AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1122      if (mgr) mgr->AddStatisticsMsg(msg);
1123   }
1124 }
1125
1126 Long64_t AliPhysicsSelection::Merge(TCollection* list)
1127 {
1128   // Merge a list of AliMultiplicityCorrection objects with this (needed for
1129   // PROOF).
1130   // Returns the number of merged objects (including this).
1131
1132   if (!list)
1133     return 0;
1134
1135   if (list->IsEmpty())
1136     return 1;
1137
1138   TIterator* iter = list->MakeIterator();
1139   TObject* obj;
1140
1141   // collections of all histograms
1142   const Int_t nHists = 9;
1143   TList collections[nHists];
1144
1145   Int_t count = 0;
1146   while ((obj = iter->Next())) {
1147
1148     AliPhysicsSelection* entry = dynamic_cast<AliPhysicsSelection*> (obj);
1149     if (entry == 0)
1150       continue;
1151     // Update run number. If this one is not initialized (-1) take the one from
1152     // the next physics selection to be merged with. In case of 2 different run
1153     // numbers issue a warning (should physics selections from different runs be
1154     // merged together) A.G.
1155     Int_t currentRun = entry->GetCurrentRun();
1156     // Nothing to merge with since run number was not initialized.
1157     if (currentRun < 0) continue;
1158     if (fCurrentRun < 0)
1159     {
1160       fCurrentRun = currentRun;
1161       fMC = entry->fMC;
1162       for (Int_t i=0; i<entry->fCollTrigClasses.GetEntries(); i++)
1163         fCollTrigClasses.Add(entry->fCollTrigClasses.At(i)->Clone());
1164       for (Int_t i=0; i<entry->fBGTrigClasses.GetEntries(); i++)
1165         fBGTrigClasses.Add(entry->fBGTrigClasses.At(i)->Clone());
1166     }
1167     if (fCurrentRun != currentRun)
1168        AliWarning(Form("Current run %d not matching the one to be merged with %d", fCurrentRun, currentRun));
1169
1170     if (entry->fTriggerAnalysis.GetEntries() > 0)
1171       collections[0].Add(&(entry->fTriggerAnalysis));
1172     if (entry->fHistStatistics[0])
1173       collections[1].Add(entry->fHistStatistics[0]);
1174     if (entry->fHistStatistics[1])
1175       collections[2].Add(entry->fHistStatistics[1]);
1176     if (entry->fHistBunchCrossing)
1177       collections[3].Add(entry->fHistBunchCrossing);
1178     if (entry->fHistTriggerPattern)
1179       collections[4].Add(entry->fHistTriggerPattern);
1180     if (entry->fBackgroundIdentification)
1181       collections[5].Add(entry->fBackgroundIdentification);
1182
1183     count++;
1184   }
1185
1186   if (fTriggerAnalysis.GetEntries() == 0 && collections[0].GetEntries() > 0)
1187   {
1188     TList* firstList = (TList*) collections[0].First();
1189     for (Int_t i=0; i<firstList->GetEntries(); i++)
1190       fTriggerAnalysis.Add(firstList->At(i)->Clone());
1191
1192     collections[0].RemoveAt(0);
1193   }
1194   fTriggerAnalysis.Merge(&collections[0]);
1195
1196   // if this instance is empty (not initialized) nothing would be merged here --> copy first entry
1197   if (!fHistStatistics[0] && collections[1].GetEntries() > 0)
1198   {
1199     fHistStatistics[0] = (TH2F*) collections[1].First()->Clone();
1200     collections[1].RemoveAt(0);
1201   }
1202   if (fHistStatistics[0])
1203     fHistStatistics[0]->Merge(&collections[1]);
1204
1205   if (!fHistStatistics[1] && collections[2].GetEntries() > 0)
1206   {
1207     fHistStatistics[1] = (TH2F*) collections[2].First()->Clone();
1208     collections[2].RemoveAt(0);
1209   }
1210   if (fHistStatistics[1])
1211     fHistStatistics[1]->Merge(&collections[2]);
1212
1213   if (!fHistBunchCrossing && collections[3].GetEntries() > 0)
1214   {
1215     fHistBunchCrossing = (TH2F*) collections[3].First()->Clone();
1216     collections[3].RemoveAt(0);
1217   }
1218   if (fHistBunchCrossing)
1219     fHistBunchCrossing->Merge(&collections[3]);
1220
1221   if (!fHistTriggerPattern && collections[4].GetEntries() > 0)
1222   {
1223     fHistTriggerPattern = (TH1F*) collections[4].First()->Clone();
1224     collections[4].RemoveAt(0);
1225   }
1226   if (fHistTriggerPattern)
1227     fHistTriggerPattern->Merge(&collections[4]);
1228
1229   if (!fBackgroundIdentification && collections[5].GetEntries() > 0)
1230   {
1231     fBackgroundIdentification = (AliAnalysisCuts*) collections[5].First()->Clone();
1232     collections[5].RemoveAt(0);
1233   }
1234   if (fBackgroundIdentification)
1235     fBackgroundIdentification->Merge(&collections[5]);
1236
1237   delete iter;
1238
1239   return count+1;
1240 }
1241
1242 void AliPhysicsSelection::SaveHistograms(const char* folder) const
1243 {
1244   // write histograms to current directory
1245
1246   if (!fHistStatistics[0] || !fHistStatistics[1])
1247     return;
1248
1249   if (folder)
1250   {
1251     gDirectory->mkdir(folder);
1252     gDirectory->cd(folder);
1253   }
1254
1255
1256   // Fill the last rows of fHistStatistics before saving
1257   if (fComputeBG) {
1258     Int_t triggerScheme = GetTriggerScheme(UInt_t(fCurrentRun));
1259     if(triggerScheme != 1){
1260       AliWarning("BG estimate only supported for trigger scheme \"1\" (CINT1 suite)");
1261     } else if (fUseMuonTriggers) {
1262       AliWarning("BG estimate with muon triggers to be implemented");
1263     } else {
1264       AliInfo("BG estimate assumes that for a given run you only have A and C triggers separately or"
1265               " toghether as a AC class! Make sure this assumption holds in your case");
1266
1267       // use an anum for the different trigger classes, to make loops easier to read
1268       enum {kClassB =0 , kClassA, kClassC, kClassAC, kClassE, kNClasses};
1269       const char * classFlags[] = {"B", "A", "C", "AC", "E"}; // labels
1270
1271       UInt_t * rows[kNClasses] = {0}; // Array of matching rows
1272       Int_t nrows[kNClasses] = {0};
1273       // Get rows matching the requested trigger bits for all trigger classes
1274       for(Int_t iTrigClass = 0; iTrigClass < kNClasses; iTrigClass++){
1275         nrows[iTrigClass] = GetStatRow(classFlags[iTrigClass],fComputeBG,&rows[iTrigClass]);
1276       }
1277
1278       // 0. Determine the ratios of triggers E/B, A/B, C/B from the stat histogram
1279       // Those are used to rescale the different classes to the same number of bx ids
1280       // TODO: pass names of the rows for B, CA and E and look names of the rows. How do I handle the case in which both AC are in the same row?
1281       Int_t nBXIds[kNClasses] = {0};
1282       cout <<"Computing BG:" << endl;
1283
1284       for(Int_t iTrigClass = 0; iTrigClass < kNClasses; iTrigClass++){
1285         for(Int_t irow = 0; irow < nrows[iTrigClass]; irow++) {
1286           if(irow==0) cout << "- Class " << classFlags[iTrigClass] << endl;
1287           for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++) {
1288             if (fHistBunchCrossing->GetBinContent(j, rows[iTrigClass][irow]) > 0) {
1289               nBXIds[iTrigClass]++;
1290             }
1291           }
1292           if(nBXIds[iTrigClass]>0) cout << "   Using row " << rows[iTrigClass][irow] <<  ": "
1293                                         << fHistBunchCrossing->GetYaxis()->GetBinLabel(rows[iTrigClass][irow])
1294                                         << " (nBXID "<< nBXIds[iTrigClass] << ")"<< endl;
1295
1296         }
1297
1298       }
1299
1300       Float_t ratioToB[kNClasses];
1301       ratioToB[kClassE]  = nBXIds[kClassE]  >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassE]   : 0;
1302       ratioToB[kClassA]  = nBXIds[kClassA]  >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassA]   : 0;
1303       ratioToB[kClassC]  = nBXIds[kClassC]  >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassC]   : 0;
1304       ratioToB[kClassAC] = nBXIds[kClassAC] >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassAC]  : 0;
1305       Printf("Ratio between the BX ids in the different trigger classes:");
1306       Printf("  B/E  = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassE], ratioToB[kClassE] );
1307       Printf("  B/A  = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassA], ratioToB[kClassA] );
1308       Printf("  B/C  = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassC], ratioToB[kClassC] );
1309       Printf("  B/AC = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassAC],ratioToB[kClassAC]);
1310       Int_t nHistStat = 2;
1311
1312       // 1. loop over all cols
1313       for(Int_t iHistStat = 0; iHistStat < nHistStat; iHistStat++){
1314         Int_t ncol = fHistStatistics[iHistStat]->GetNbinsX();
1315         Float_t good1 = 0;
1316         for(Int_t icol = 1; icol <= ncol; icol++) {
1317           Int_t nEvents[kNClasses] = {0}; // number of events should be reset at every column
1318           // For all trigger classes, add up over row matching trigger mask (as selected before)
1319           for(Int_t iTrigClass = 0; iTrigClass < kNClasses; iTrigClass++){
1320             for(Int_t irow = 0; irow < nrows[iTrigClass]; irow++) {
1321               nEvents[iTrigClass] += (Int_t) fHistStatistics[iHistStat]->GetBinContent(icol,rows[iTrigClass][irow]);
1322             }
1323             //      cout << "Events " << classFlags[iTrigClass] << " ("<<icol<<") " << nEvents[iTrigClass] << endl;
1324           }
1325           if (nEvents[kClassB]>0) {
1326             Float_t acc  = ratioToB[kClassE]*nEvents[kClassE];
1327             Double_t acc_err = TMath::Sqrt(ratioToB[kClassE]*ratioToB[kClassE]*nEvents[kClassE]);
1328             //      Int_t bg   = cint1A + cint1C - 2*acc;
1329
1330              // Assuming that for a given class the triggers are either recorded as A+C or AC
1331             Float_t bg   = nEvents[kClassAC] > 0 ?
1332               fBIFactorAC*(ratioToB[kClassAC]*nEvents[kClassAC] - 2*acc):
1333               fBIFactorA* (ratioToB[kClassA]*nEvents[kClassA]-acc) +
1334               fBIFactorC* (ratioToB[kClassC]*nEvents[kClassC]-acc) ;
1335
1336             // cout << "-----------------------" << endl;
1337             // cout << "Factors: " << fBIFactorA << " " << fBIFactorC << " " << fBIFactorAC << endl;
1338             // cout << "Ratios: "  << ratioToB[kClassA] << " " << ratioToB[kClassC] << " " << ratioToB[kClassAC] << endl;
1339             // cout << "Evts:   "  << nEvents[kClassA] << " " << nEvents[kClassC] << " " << nEvents[kClassAC] << " " <<  nEvents[kClassB] << endl;
1340             // cout << "Acc: " << acc << endl;
1341             // cout << "BG: " << bg << endl;
1342             // cout  << "  " <<   fBIFactorA* (ratioToB[kClassA]*nEvents[kClassA]-acc) <<endl;
1343             // cout  << "  " <<   fBIFactorC* (ratioToB[kClassC]*nEvents[kClassC]-acc) <<endl;
1344             // cout  << "  " <<   fBIFactorAC*(ratioToB[kClassAC]*nEvents[kClassAC] - 2*acc) << endl;
1345             // cout << "-----------------------" << endl;
1346
1347             Float_t good = Float_t(nEvents[kClassB]) - bg - acc;
1348             if (icol ==1) good1 = good;
1349             //      Float_t errGood     = TMath::Sqrt(2*(nEvents[kClassA]+nEvents[kClassC]+nEvents[kClassE]));// Error on the number of goods assuming only bg fluctuates
1350             //      DeltaG^2 = B + FA^2 A + FC^2 C + Ratio^2 (FA+FC-1)^2 E.
1351             Float_t errGood     = nEvents[kClassAC] > 0 ?
1352               TMath::Sqrt( nEvents[kClassB] +
1353                            fBIFactorAC*fBIFactorAC*ratioToB[kClassAC]*ratioToB[kClassAC]*nEvents[kClassAC] +
1354                            ratioToB[kClassE] * ratioToB[kClassE] *
1355                            (fBIFactorAC - 1)*(fBIFactorAC - 1)*nEvents[kClassE]) :
1356               TMath::Sqrt( nEvents[kClassB] +
1357                            fBIFactorA*fBIFactorA*ratioToB[kClassA]*ratioToB[kClassA]*nEvents[kClassA] +
1358                            fBIFactorC*fBIFactorC*ratioToB[kClassC]*ratioToB[kClassC]*nEvents[kClassC] +
1359                            ratioToB[kClassE] * ratioToB[kClassE] *
1360                            (fBIFactorA + fBIFactorC - 1)*(fBIFactorA + fBIFactorC - 1)*nEvents[kClassE]);
1361
1362             Float_t errBG = nEvents[kClassAC] > 0 ?
1363               TMath::Sqrt(fBIFactorAC*fBIFactorAC*ratioToB[kClassAC]*ratioToB[kClassAC]*nEvents[kClassAC]+
1364                           4*ratioToB[kClassE]*ratioToB[kClassE]*(fBIFactorAC*fBIFactorAC)*nEvents[kClassE]) :
1365               TMath::Sqrt(fBIFactorA*fBIFactorA*ratioToB[kClassA]*ratioToB[kClassA]*nEvents[kClassA]+
1366                           fBIFactorC*fBIFactorC*ratioToB[kClassC]*ratioToB[kClassC]*nEvents[kClassC]+
1367                           ratioToB[kClassE]*ratioToB[kClassE]*(fBIFactorA+fBIFactorC)*(fBIFactorA+fBIFactorC)*nEvents[kClassE]);
1368
1369
1370             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowBG,bg);
1371             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowBG,errBG);
1372             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowAcc,acc);
1373             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowAcc,acc_err);
1374             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowGood,good);
1375             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowGood,errGood);
1376
1377 #ifdef VERBOSE_STAT
1378             //kStatRowBG=0,kStatRowAcc,kStatRowBGFrac,kStatRowAccFrac,kStatRowErrGoodFrac,kStatRowGoodFrac,kStatRowGood,kStatRowErrGood
1379             Float_t accFrac   = Float_t(acc) / nEvents[kClassB]  *100;
1380             Float_t errAccFrac= Float_t(acc_err) / nEvents[kClassB]  *100;
1381             Float_t bgFrac    = Float_t(bg)  / nEvents[kClassB]  *100;
1382             Float_t goodFrac  = Float_t(good)  / good1 *100;
1383             Float_t errGoodFrac = errGood/good1 * 100;
1384             Float_t errFracBG = bg > 0 ? TMath::Sqrt((errBG/bg)*(errBG/bg) + 1/nEvents[kClassB])*bgFrac : 0;
1385             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowBGFrac,bgFrac);
1386             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowBGFrac,errFracBG);
1387             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowAccFrac,accFrac);
1388             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowAccFrac,errAccFrac);
1389             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowGoodFrac,goodFrac);
1390             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowErrGoodFrac,errGoodFrac);
1391             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowErrGood,errGood);
1392 #endif
1393           }
1394         }
1395       }
1396       for (Int_t iTrigClass = 0; iTrigClass < kNClasses; iTrigClass++){
1397         delete [] rows[iTrigClass];
1398       }
1399     }
1400   }
1401
1402   fHistStatistics[0]->Write();
1403   fHistStatistics[1]->Write();
1404   if(fHistBunchCrossing ) fHistBunchCrossing ->Write();
1405   if(fHistTriggerPattern) fHistTriggerPattern->Write();
1406
1407   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
1408   for (Int_t i=0; i < count; i++)
1409   {
1410     TString triggerClass = "trigger_histograms_";
1411     if (i < fCollTrigClasses.GetEntries())
1412       triggerClass += ((TObjString*) fCollTrigClasses.At(i))->String();
1413     else
1414       triggerClass += ((TObjString*) fBGTrigClasses.At(i - fCollTrigClasses.GetEntries()))->String();
1415
1416     gDirectory->mkdir(triggerClass);
1417     gDirectory->cd(triggerClass);
1418
1419     static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i))->SaveHistograms();
1420
1421     gDirectory->cd("..");
1422   }
1423
1424   if (fBackgroundIdentification)
1425   {
1426     gDirectory->mkdir("background_identification");
1427     gDirectory->cd("background_identification");
1428
1429     fBackgroundIdentification->GetOutput()->Write();
1430
1431     gDirectory->cd("..");
1432   }
1433
1434   if (folder)
1435     gDirectory->cd("..");
1436 }
1437
1438 Int_t AliPhysicsSelection::GetStatRow(const char * triggerBXClass, UInt_t offlineTriggerType, UInt_t ** rowIDs) const {
1439   // Puts inside the array rowIDs the row number for a given offline
1440   // trigger in a given bx class. Returns the total number of lines
1441   // matching the selection
1442   // triggerBXClass can be either "A", "AC", "B" or "E"
1443   // offlineTriggerType is one of the types defined in AliVEvent
1444   // User should delete rowIDs if no longer needed
1445
1446   if(!fHistStatistics[0]) {
1447     AliWarning("Not initialized, returning 0");
1448     return 0;
1449   }
1450   const Int_t nrows = fHistStatistics[0]->GetNbinsY();
1451
1452   // allocate memory for at maximum nrows
1453   Int_t nMatches = 0;
1454   (*rowIDs) = new UInt_t[nrows];
1455
1456   // Build regular expression. look for a +, followed by the beginning
1457   // of a word. Within the word, look for the class id after any
1458   // number of any char, but at most one dash ("[^-]*-?"), followed by
1459   // a - and then any char (".*") and at the class id ("\\&(\\d)")
1460   // The class id is stored.
1461   // WARNING: please check this if the trigger classes change
1462   TPRegexp re(Form("\\+\\b[^-]*-?%s-.*\\&(\\d)",triggerBXClass));
1463   // Loop over rows and find matching ones:
1464   for(Int_t irow = 1; irow <= nrows; irow++){
1465     TObjArray * matches = re.MatchS(fHistStatistics[0]->GetYaxis()->GetBinLabel(irow));
1466     if (matches->GetEntries()) {
1467       TString s = ((TObjString*)matches->At(1))->GetString();
1468       if(UInt_t(s.Atoi()) & offlineTriggerType) { // bitwise comparison with the requested mask
1469         //      cout << "Marching " << s.Data() << " " << offlineTriggerType << " " << fHistStatistics[0]->GetYaxis()->GetBinLabel(irow) << endl;
1470         (*rowIDs)[nMatches] = irow;
1471         nMatches++;
1472       }
1473     }
1474     delete matches;
1475   }
1476
1477   return nMatches;
1478 }
1479
1480 void AliPhysicsSelection::SetBIFactors(const AliESDEvent * aESD) {
1481   // Set factors for realtive bunch intesities
1482   if(!aESD) {
1483     AliFatal("ESD not given");
1484   }
1485   Int_t run = aESD->GetRunNumber();
1486   if (run > 105268) {
1487     // intensities stored in the ESDs
1488     const AliESDRun* esdRun = aESD->GetESDRun();
1489     Double_t intAB = esdRun->GetMeanIntensityIntecting(0);
1490     Double_t intCB = esdRun->GetMeanIntensityIntecting(1);
1491     Double_t intAA = esdRun->GetMeanIntensityNonIntecting(0);
1492     Double_t intCC = esdRun->GetMeanIntensityNonIntecting(1);
1493
1494     // cout << "INT " <<intAB <<endl;
1495     // cout << "INT " <<intCB <<endl;
1496     // cout << "INT " <<intAA <<endl;
1497     // cout << "INT " <<intCC <<endl;
1498
1499     if (intAB > -1 && intAA > -1) {
1500       fBIFactorA = intAB/intAA;
1501     } else {
1502       AliWarning("Cannot set fBIFactorA, assuming 1");
1503     }
1504
1505     if (intCB > -1 && intCC > -1) {
1506       fBIFactorC = intCB/intCC;
1507     } else {
1508       AliWarning("Cannot set fBIFactorC, assuming 1");
1509     }
1510
1511     if (intAB > -1 && intAA > -1 &&
1512         intCB > -1 && intCC > -1) {
1513       fBIFactorAC = (intAB+intCB)/(intAA+intCC);
1514     } else {
1515       AliWarning("Cannot set fBIFactorAC, assuming 1");
1516     }
1517
1518   }
1519   else {
1520     // First runs. Intensities hardcoded
1521     switch(run) {
1522     case 104155:
1523       fBIFactorA = 0.961912722908;
1524       fBIFactorC = 1.04992336081;
1525       break;
1526     case 104157:
1527       fBIFactorA = 0.947312854998;
1528       fBIFactorC = 1.01599706417;
1529       break;
1530     case 104159:
1531       fBIFactorA = 0.93659320151;
1532       fBIFactorC = 0.98580804207;
1533       break;
1534     case 104160:
1535       fBIFactorA = 0.929664189926;
1536       fBIFactorC = 0.963467679851;
1537       break;
1538     case 104315:
1539       fBIFactorA = 1.08939104979;
1540       fBIFactorC = 0.931113921925;
1541       break;
1542     case 104316:
1543       fBIFactorA = 1.08351880974;
1544       fBIFactorC = 0.916068345845;
1545       break;
1546     case 104320:
1547       fBIFactorA = 1.07669281245;
1548       fBIFactorC = 0.876818744763;
1549       break;
1550     case 104321:
1551       fBIFactorA = 1.00971079602;
1552       fBIFactorC = 0.773781299076;
1553       break;
1554     case 104792:
1555       fBIFactorA = 0.787215863962;
1556       fBIFactorC = 0.778253173071;
1557       break;
1558     case 104793:
1559       fBIFactorA = 0.692211363661;
1560       fBIFactorC = 0.733152456667;
1561       break;
1562     case 104799:
1563       fBIFactorA = 1.04027825161;
1564       fBIFactorC = 1.00530825942;
1565       break;
1566     case 104800:
1567       fBIFactorA = 1.05309910671;
1568       fBIFactorC = 1.00376801855;
1569       break;
1570     case 104801:
1571       fBIFactorA = 1.0531231922;
1572       fBIFactorC = 0.992439666758;
1573       break;
1574     case 104802:
1575       fBIFactorA = 1.04191478134;
1576       fBIFactorC = 0.979368585208;
1577       break;
1578     case 104803:
1579       fBIFactorA = 1.03121314094;
1580       fBIFactorC = 0.973379962609;
1581       break;
1582     case 104824:
1583       fBIFactorA = 0.969945926722;
1584       fBIFactorC = 0.39549745806;
1585       break;
1586     case 104825:
1587       fBIFactorA = 0.968627213937;
1588       fBIFactorC = 0.310100412205;
1589       break;
1590     case 104841:
1591       fBIFactorA = 0.991601393212;
1592       fBIFactorC = 0.83762204722;
1593       break;
1594     case 104845:
1595       fBIFactorA = 0.98040863886;
1596       fBIFactorC = 0.694824205793;
1597       break;
1598     case 104867:
1599       fBIFactorA = 1.10646173412;
1600       fBIFactorC = 0.841407246916;
1601       break;
1602     case 104876:
1603       fBIFactorA = 1.12063452421;
1604       fBIFactorC = 0.78726542895;
1605       break;
1606     case 104890:
1607       fBIFactorA = 1.02346137453;
1608       fBIFactorC = 1.03355663595;
1609       break;
1610     case 104892:
1611       fBIFactorA = 1.05406025913;
1612       fBIFactorC = 1.00029166135;
1613       break;
1614     case 105143:
1615       fBIFactorA = 0.947343384349;
1616       fBIFactorC = 0.972637444408;
1617       break;
1618     case 105160:
1619       fBIFactorA = 0.908854622177;
1620       fBIFactorC = 0.958851103977;
1621       break;
1622     case 105256:
1623       fBIFactorA = 0.810076150206;
1624       fBIFactorC = 0.884663561883;
1625       break;
1626     case 105257:
1627       fBIFactorA = 0.80974912303;
1628       fBIFactorC = 0.878859123479;
1629       break;
1630     case 105268:
1631       fBIFactorA = 0.809052110679;
1632       fBIFactorC = 0.87233890989;
1633       break;
1634     default:
1635       fBIFactorA = 1;
1636       fBIFactorC = 1;
1637     }
1638   }
1639
1640 }