ANALYSIS/AliPhysicsSelection.cxx

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