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