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