ANALYSIS/AliPhysicsSelection.cxx

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