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