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   TString msg;
 934   Printf("Configuration initialized for run %d (MC: %d):", fCurrentRun, fMC);
 935   msg += Form("Configuration initialized for run %d (MC: %d):\n", fCurrentRun, fMC);
 936
 937   Printf("Collision trigger classes:");
 938   for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
 939     Printf("%s", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
 940
 941   Printf("Background trigger classes:");
 942   for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
 943     Printf("%s", ((TObjString*) fBGTrigClasses.At(i))->String().Data());
 944
 945   AliTriggerAnalysis* triggerAnalysis = dynamic_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(0));
 946
 947   if (triggerAnalysis)
 948   {
 949     if (triggerAnalysis->GetV0TimeOffset() > 0)
 950       Printf("V0 time offset active: %.2f ns", triggerAnalysis->GetV0TimeOffset());
 951
 952     Printf("\nTotal available events:");
 953
 954     triggerAnalysis->PrintTriggerClasses();
 955   }
 956
 957   if (fHistStatistics[kStatIdxAll])
 958   {
 959     for (Int_t i=0; i<fCollTrigClasses.GetEntries(); i++)
 960     {
 961       Printf("\nSelection statistics for collision trigger %s:", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
 962       msg += Form("\nSelection statistics for collision trigger %s:\n", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
 963
 964       Printf("Total events with correct trigger class: %d", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(1, i+1));
 965       msg += Form("Total events with correct trigger class: %d\n", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(1, i+1));
 966       Printf("Selected collision candidates: %d", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(fHistStatistics[kStatIdxAll]->GetXaxis()->FindBin("Accepted"), i+1));
 967       msg += Form("Selected collision candidates: %d\n", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(fHistStatistics[kStatIdxAll]->GetXaxis()->FindBin("Accepted"), i+1));
 968     }
 969   }
 970
 971   if (fHistBunchCrossing)
 972   {
 973     Printf("\nBunch crossing statistics:");
 974     msg += "\nBunch crossing statistics:\n";
 975
 976     for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
 977     {
 978       TString str;
 979       str.Form("Trigger %s has accepted events in the bunch crossings: ", fHistBunchCrossing->GetYaxis()->GetBinLabel(i));
 980
 981       for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
 982         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
 983           str += Form("%d, ", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
 984
 985       Printf("%s", str.Data());
 986       msg += str;
 987       msg += "\n";
 988     }
 989
 990     for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
 991     {
 992       Int_t countColl = 0;
 993       Int_t countBG = 0;
 994       for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
 995       {
 996         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
 997         {
 998           if (fCollTrigClasses.FindObject(fHistBunchCrossing->GetYaxis()->GetBinLabel(i)))
 999             countColl++;
1000           if (fBGTrigClasses.FindObject(fHistBunchCrossing->GetYaxis()->GetBinLabel(i)))
1001             countBG++;
1002         }
1003       }
1004       if (countColl > 0 && countBG > 0)
1005         Printf("WARNING: Bunch crossing %d has collision and BG trigger classes active. Check BPTX functioning for this run!", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
1006     }
1007   }
1008
1009   if (fUsingCustomClasses)
1010     Printf("WARNING: Using custom trigger classes!");
1011   if (fSkipTriggerClassSelection)
1012     Printf("WARNING: Skipping trigger class selection!");
1013   if (fSkipV0)
1014     Printf("WARNING: Ignoring V0 information in selection");
1015   if(!fBin0CallBack)
1016     Printf("WARNING: Callback not set: will not fill the statistics for the bin 0");
1017   TString opt(option);
1018   opt.ToUpper();
1019   if (opt == "STAT") {
1020      AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1021      if (mgr) mgr->AddStatisticsMsg(msg);
1022   }
1023 }
1024
1025 Long64_t AliPhysicsSelection::Merge(TCollection* list)
1026 {
1027   // Merge a list of AliMultiplicityCorrection objects with this (needed for
1028   // PROOF).
1029   // Returns the number of merged objects (including this).
1030
1031   if (!list)
1032     return 0;
1033
1034   if (list->IsEmpty())
1035     return 1;
1036
1037   TIterator* iter = list->MakeIterator();
1038   TObject* obj;
1039
1040   // collections of all histograms
1041   const Int_t nHists = 9;
1042   TList collections[nHists];
1043
1044   Int_t count = 0;
1045   while ((obj = iter->Next())) {
1046
1047     AliPhysicsSelection* entry = dynamic_cast<AliPhysicsSelection*> (obj);
1048     if (entry == 0)
1049       continue;
1050     // Update run number. If this one is not initialized (-1) take the one from
1051     // the next physics selection to be merged with. In case of 2 different run
1052     // numbers issue a warning (should physics selections from different runs be
1053     // merged together) A.G.
1054     Int_t currentRun = entry->GetCurrentRun();
1055     // Nothing to merge with since run number was not initialized.
1056     if (currentRun < 0) continue;
1057     if (fCurrentRun < 0) fCurrentRun = currentRun;
1058     if (fCurrentRun != currentRun)
1059        AliWarning(Form("Current run %d not matching the one to be merged with %d", fCurrentRun, currentRun));
1060
1061     collections[0].Add(&(entry->fTriggerAnalysis));
1062     if (entry->fHistStatistics[0])
1063       collections[1].Add(entry->fHistStatistics[0]);
1064     if (entry->fHistStatistics[1])
1065       collections[2].Add(entry->fHistStatistics[1]);
1066     if (entry->fHistBunchCrossing)
1067       collections[3].Add(entry->fHistBunchCrossing);
1068     if (entry->fHistTriggerPattern)
1069       collections[4].Add(entry->fHistTriggerPattern);
1070     if (entry->fBackgroundIdentification)
1071       collections[5].Add(entry->fBackgroundIdentification);
1072
1073     count++;
1074   }
1075
1076   fTriggerAnalysis.Merge(&collections[0]);
1077   if (fHistStatistics[0])
1078     fHistStatistics[0]->Merge(&collections[1]);
1079   if (fHistStatistics[1])
1080     fHistStatistics[1]->Merge(&collections[2]);
1081   if (fHistBunchCrossing)
1082     fHistBunchCrossing->Merge(&collections[3]);
1083   if (fHistTriggerPattern)
1084     fHistTriggerPattern->Merge(&collections[4]);
1085   if (fBackgroundIdentification)
1086     fBackgroundIdentification->Merge(&collections[5]);
1087
1088   delete iter;
1089
1090   return count+1;
1091 }
1092
1093 void AliPhysicsSelection::SaveHistograms(const char* folder) const
1094 {
1095   // write histograms to current directory
1096
1097   if (!fHistStatistics[0] || !fHistStatistics[1])
1098     return;
1099
1100   if (folder)
1101   {
1102     gDirectory->mkdir(folder);
1103     gDirectory->cd(folder);
1104   }
1105
1106
1107   // Fill the last rows of fHistStatistics before saving
1108   if (fComputeBG) {
1109     Int_t triggerScheme = GetTriggerScheme(UInt_t(fCurrentRun));
1110     if(triggerScheme != 1){
1111       AliWarning("BG estimate only supported for trigger scheme \"1\" (CINT1 suite)");
1112     } else if (fUseMuonTriggers) {
1113       AliWarning("BG estimate with muon triggers to be implemented");
1114     } else {
1115       AliInfo("BG estimate assumes that for a given run you only have A and C triggers separately or"
1116               " toghether as a AC class! Make sure this assumption holds in your case");
1117
1118       // use an anum for the different trigger classes, to make loops easier to read
1119       enum {kClassB =0 , kClassA, kClassC, kClassAC, kClassE, kNClasses};
1120       const char * classFlags[] = {"B", "A", "C", "AC", "E"}; // labels
1121
1122       UInt_t * rows[kNClasses] = {0}; // Array of matching rows
1123       Int_t nrows[kNClasses] = {0};
1124       // Get rows matching the requested trigger bits for all trigger classes
1125       for(Int_t iTrigClass = 0; iTrigClass < kNClasses; iTrigClass++){
1126         nrows[iTrigClass] = GetStatRow(classFlags[iTrigClass],fComputeBG,&rows[iTrigClass]);
1127       }
1128
1129       // 0. Determine the ratios of triggers E/B, A/B, C/B from the stat histogram
1130       // Those are used to rescale the different classes to the same number of bx ids
1131       // 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?
1132       Int_t nBXIds[kNClasses] = {0};
1133       cout <<"Computing BG:" << endl;
1134
1135       for(Int_t iTrigClass = 0; iTrigClass < kNClasses; iTrigClass++){
1136         for(Int_t irow = 0; irow < nrows[iTrigClass]; irow++) {
1137           if(irow==0) cout << "- Class " << classFlags[iTrigClass] << endl;
1138           for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++) {
1139             if (fHistBunchCrossing->GetBinContent(j, rows[iTrigClass][irow]) > 0) {
1140               nBXIds[iTrigClass]++;
1141             }
1142           }
1143           if(nBXIds[iTrigClass]>0) cout << "   Using row " << rows[iTrigClass][irow] <<  ": "
1144                                         << fHistBunchCrossing->GetYaxis()->GetBinLabel(rows[iTrigClass][irow])
1145                                         << " (nBXID "<< nBXIds[iTrigClass] << ")"<< endl;
1146
1147         }
1148
1149       }
1150
1151       Float_t ratioToB[kNClasses];
1152       ratioToB[kClassE]  = nBXIds[kClassE]  >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassE]   : 0;
1153       ratioToB[kClassA]  = nBXIds[kClassA]  >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassA]   : 0;
1154       ratioToB[kClassC]  = nBXIds[kClassC]  >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassC]   : 0;
1155       ratioToB[kClassAC] = nBXIds[kClassAC] >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassAC]  : 0;
1156       Printf("Ratio between the BX ids in the different trigger classes:");
1157       Printf("  B/E  = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassE], ratioToB[kClassE] );
1158       Printf("  B/A  = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassA], ratioToB[kClassA] );
1159       Printf("  B/C  = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassC], ratioToB[kClassC] );
1160       Printf("  B/AC = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassAC],ratioToB[kClassAC]);
1161       Int_t nHistStat = 2;
1162
1163       // 1. loop over all cols
1164       for(Int_t iHistStat = 0; iHistStat < nHistStat; iHistStat++){
1165         Int_t ncol = fHistStatistics[iHistStat]->GetNbinsX();
1166         Float_t good1 = 0;
1167         for(Int_t icol = 1; icol <= ncol; icol++) {
1168           Int_t nEvents[kNClasses] = {0}; // number of events should be reset at every column
1169           // For all trigger classes, add up over row matching trigger mask (as selected before)
1170           for(Int_t iTrigClass = 0; iTrigClass < kNClasses; iTrigClass++){
1171             for(Int_t irow = 0; irow < nrows[iTrigClass]; irow++) {
1172               nEvents[iTrigClass] += (Int_t) fHistStatistics[iHistStat]->GetBinContent(icol,rows[iTrigClass][irow]);
1173             }
1174             //      cout << "Events " << classFlags[iTrigClass] << " ("<<icol<<") " << nEvents[iTrigClass] << endl;
1175           }
1176           if (nEvents[kClassB]>0) {
1177             Float_t acc  = ratioToB[kClassE]*nEvents[kClassE];
1178             Double_t acc_err = TMath::Sqrt(ratioToB[kClassE]*ratioToB[kClassE]*nEvents[kClassE]);
1179             //      Int_t bg   = cint1A + cint1C - 2*acc;
1180
1181              // Assuming that for a given class the triggers are either recorded as A+C or AC
1182             Float_t bg   = nEvents[kClassAC] > 0 ?
1183               fBIFactorAC*(ratioToB[kClassAC]*nEvents[kClassAC] - 2*acc):
1184               fBIFactorA* (ratioToB[kClassA]*nEvents[kClassA]-acc) +
1185               fBIFactorC* (ratioToB[kClassC]*nEvents[kClassC]-acc) ;
1186
1187             // cout << "-----------------------" << endl;
1188             // cout << "Factors: " << fBIFactorA << " " << fBIFactorC << " " << fBIFactorAC << endl;
1189             // cout << "Ratios: "  << ratioToB[kClassA] << " " << ratioToB[kClassC] << " " << ratioToB[kClassAC] << endl;
1190             // cout << "Evts:   "  << nEvents[kClassA] << " " << nEvents[kClassC] << " " << nEvents[kClassAC] << " " <<  nEvents[kClassB] << endl;
1191             // cout << "Acc: " << acc << endl;
1192             // cout << "BG: " << bg << endl;
1193             // cout  << "  " <<   fBIFactorA* (ratioToB[kClassA]*nEvents[kClassA]-acc) <<endl;
1194             // cout  << "  " <<   fBIFactorC* (ratioToB[kClassC]*nEvents[kClassC]-acc) <<endl;
1195             // cout  << "  " <<   fBIFactorAC*(ratioToB[kClassAC]*nEvents[kClassAC] - 2*acc) << endl;
1196             // cout << "-----------------------" << endl;
1197
1198             Float_t good = Float_t(nEvents[kClassB]) - bg - acc;
1199             if (icol ==1) good1 = good;
1200             //      Float_t errGood     = TMath::Sqrt(2*(nEvents[kClassA]+nEvents[kClassC]+nEvents[kClassE]));// Error on the number of goods assuming only bg fluctuates
1201             //      DeltaG^2 = B + FA^2 A + FC^2 C + Ratio^2 (FA+FC-1)^2 E.
1202             Float_t errGood     = nEvents[kClassAC] > 0 ?
1203               TMath::Sqrt( nEvents[kClassB] +
1204                            fBIFactorAC*fBIFactorAC*ratioToB[kClassAC]*ratioToB[kClassAC]*nEvents[kClassAC] +
1205                            ratioToB[kClassE] * ratioToB[kClassE] *
1206                            (fBIFactorAC - 1)*(fBIFactorAC - 1)*nEvents[kClassE]) :
1207               TMath::Sqrt( nEvents[kClassB] +
1208                            fBIFactorA*fBIFactorA*ratioToB[kClassA]*ratioToB[kClassA]*nEvents[kClassA] +
1209                            fBIFactorC*fBIFactorC*ratioToB[kClassC]*ratioToB[kClassC]*nEvents[kClassC] +
1210                            ratioToB[kClassE] * ratioToB[kClassE] *
1211                            (fBIFactorA + fBIFactorC - 1)*(fBIFactorA + fBIFactorC - 1)*nEvents[kClassE]);
1212
1213             Float_t errBG = nEvents[kClassAC] > 0 ?
1214               TMath::Sqrt(fBIFactorAC*fBIFactorAC*ratioToB[kClassAC]*ratioToB[kClassAC]*nEvents[kClassAC]+
1215                           4*ratioToB[kClassE]*ratioToB[kClassE]*(fBIFactorAC*fBIFactorAC)*nEvents[kClassE]) :
1216               TMath::Sqrt(fBIFactorA*fBIFactorA*ratioToB[kClassA]*ratioToB[kClassA]*nEvents[kClassA]+
1217                           fBIFactorC*fBIFactorC*ratioToB[kClassC]*ratioToB[kClassC]*nEvents[kClassC]+
1218                           ratioToB[kClassE]*ratioToB[kClassE]*(fBIFactorA+fBIFactorC)*(fBIFactorA+fBIFactorC)*nEvents[kClassE]);
1219
1220
1221             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowBG,bg);
1222             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowBG,errBG);
1223             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowAcc,acc);
1224             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowAcc,acc_err);
1225             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowGood,good);
1226             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowGood,errGood);
1227
1228 #ifdef VERBOSE_STAT
1229             //kStatRowBG=0,kStatRowAcc,kStatRowBGFrac,kStatRowAccFrac,kStatRowErrGoodFrac,kStatRowGoodFrac,kStatRowGood,kStatRowErrGood
1230             Float_t accFrac   = Float_t(acc) / nEvents[kClassB]  *100;
1231             Float_t errAccFrac= Float_t(acc_err) / nEvents[kClassB]  *100;
1232             Float_t bgFrac    = Float_t(bg)  / nEvents[kClassB]  *100;
1233             Float_t goodFrac  = Float_t(good)  / good1 *100;
1234             Float_t errGoodFrac = errGood/good1 * 100;
1235             Float_t errFracBG = bg > 0 ? TMath::Sqrt((errBG/bg)*(errBG/bg) + 1/nEvents[kClassB])*bgFrac : 0;
1236             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowBGFrac,bgFrac);
1237             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowBGFrac,errFracBG);
1238             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowAccFrac,accFrac);
1239             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowAccFrac,errAccFrac);
1240             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowGoodFrac,goodFrac);
1241             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowErrGoodFrac,errGoodFrac);
1242             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowErrGood,errGood);
1243 #endif
1244           }
1245         }
1246       }
1247       for (Int_t iTrigClass = 0; iTrigClass < kNClasses; iTrigClass++){
1248         delete [] rows[iTrigClass];
1249       }
1250     }
1251   }
1252
1253   fHistStatistics[0]->Write();
1254   fHistStatistics[1]->Write();
1255   if(fHistBunchCrossing ) fHistBunchCrossing ->Write();
1256   if(fHistTriggerPattern) fHistTriggerPattern->Write();
1257
1258   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
1259   for (Int_t i=0; i < count; i++)
1260   {
1261     TString triggerClass = "trigger_histograms_";
1262     if (i < fCollTrigClasses.GetEntries())
1263       triggerClass += ((TObjString*) fCollTrigClasses.At(i))->String();
1264     else
1265       triggerClass += ((TObjString*) fBGTrigClasses.At(i - fCollTrigClasses.GetEntries()))->String();
1266
1267     gDirectory->mkdir(triggerClass);
1268     gDirectory->cd(triggerClass);
1269
1270     static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i))->SaveHistograms();
1271
1272     gDirectory->cd("..");
1273   }
1274
1275   if (fBackgroundIdentification)
1276   {
1277     gDirectory->mkdir("background_identification");
1278     gDirectory->cd("background_identification");
1279
1280     fBackgroundIdentification->GetOutput()->Write();
1281
1282     gDirectory->cd("..");
1283   }
1284
1285   if (folder)
1286     gDirectory->cd("..");
1287 }
1288
1289 Int_t AliPhysicsSelection::GetStatRow(const char * triggerBXClass, UInt_t offlineTriggerType, UInt_t ** rowIDs) const {
1290   // Puts inside the array rowIDs the row number for a given offline
1291   // trigger in a given bx class. Returns the total number of lines
1292   // matching the selection
1293   // triggerBXClass can be either "A", "AC", "B" or "E"
1294   // offlineTriggerType is one of the types defined in AliVEvent
1295   // User should delete rowIDs if no longer needed
1296
1297   if(!fHistStatistics[0]) {
1298     AliWarning("Not initialized, returning 0");
1299     return 0;
1300   }
1301   const Int_t nrows = fHistStatistics[0]->GetNbinsY();
1302
1303   // allocate memory for at maximum nrows
1304   Int_t nMatches = 0;
1305   (*rowIDs) = new UInt_t[nrows];
1306
1307   // Build regular expression. look for a +, followed by the beginning
1308   // of a word. Within the word, look for the class id after any
1309   // number of any char, but at most one dash ("[^-]*-?"), followed by
1310   // a - and then any char (".*") and at the class id ("\\&(\\d)")
1311   // The class id is stored.
1312   // WARNING: please check this if the trigger classes change
1313   TPRegexp re(Form("\\+\\b[^-]*-?%s-.*\\&(\\d)",triggerBXClass));
1314   // Loop over rows and find matching ones:
1315   for(Int_t irow = 1; irow <= nrows; irow++){
1316     TObjArray * matches = re.MatchS(fHistStatistics[0]->GetYaxis()->GetBinLabel(irow));
1317     if (matches->GetEntries()) {
1318       TString s = ((TObjString*)matches->At(1))->GetString();
1319       if(UInt_t(s.Atoi()) & offlineTriggerType) { // bitwise comparison with the requested mask
1320         //      cout << "Marching " << s.Data() << " " << offlineTriggerType << " " << fHistStatistics[0]->GetYaxis()->GetBinLabel(irow) << endl;
1321         (*rowIDs)[nMatches] = irow;
1322         nMatches++;
1323       }
1324     }
1325     delete matches;
1326   }
1327
1328   return nMatches;
1329 }
1330
1331 void AliPhysicsSelection::SetBIFactors(const AliESDEvent * aESD) {
1332   // Set factors for realtive bunch intesities
1333   if(!aESD) {
1334     AliFatal("ESD not given");
1335   }
1336   Int_t run = aESD->GetRunNumber();
1337   if (run > 105268) {
1338     // intensities stored in the ESDs
1339     const AliESDRun* esdRun = aESD->GetESDRun();
1340     Double_t intAB = esdRun->GetMeanIntensityIntecting(0);
1341     Double_t intCB = esdRun->GetMeanIntensityIntecting(1);
1342     Double_t intAA = esdRun->GetMeanIntensityNonIntecting(0);
1343     Double_t intCC = esdRun->GetMeanIntensityNonIntecting(1);
1344
1345     // cout << "INT " <<intAB <<endl;
1346     // cout << "INT " <<intCB <<endl;
1347     // cout << "INT " <<intAA <<endl;
1348     // cout << "INT " <<intCC <<endl;
1349
1350     if (intAB > -1 && intAA > -1) {
1351       fBIFactorA = intAB/intAA;
1352     } else {
1353       AliWarning("Cannot set fBIFactorA, assuming 1");
1354     }
1355
1356     if (intCB > -1 && intCC > -1) {
1357       fBIFactorC = intCB/intCC;
1358     } else {
1359       AliWarning("Cannot set fBIFactorC, assuming 1");
1360     }
1361
1362     if (intAB > -1 && intAA > -1 &&
1363         intCB > -1 && intCC > -1) {
1364       fBIFactorAC = (intAB+intCB)/(intAA+intCC);
1365     } else {
1366       AliWarning("Cannot set fBIFactorAC, assuming 1");
1367     }
1368
1369   }
1370   else {
1371     // First runs. Intensities hardcoded
1372     switch(run) {
1373     case 104155:
1374       fBIFactorA = 0.961912722908;
1375       fBIFactorC = 1.04992336081;
1376       break;
1377     case 104157:
1378       fBIFactorA = 0.947312854998;
1379       fBIFactorC = 1.01599706417;
1380       break;
1381     case 104159:
1382       fBIFactorA = 0.93659320151;
1383       fBIFactorC = 0.98580804207;
1384       break;
1385     case 104160:
1386       fBIFactorA = 0.929664189926;
1387       fBIFactorC = 0.963467679851;
1388       break;
1389     case 104315:
1390       fBIFactorA = 1.08939104979;
1391       fBIFactorC = 0.931113921925;
1392       break;
1393     case 104316:
1394       fBIFactorA = 1.08351880974;
1395       fBIFactorC = 0.916068345845;
1396       break;
1397     case 104320:
1398       fBIFactorA = 1.07669281245;
1399       fBIFactorC = 0.876818744763;
1400       break;
1401     case 104321:
1402       fBIFactorA = 1.00971079602;
1403       fBIFactorC = 0.773781299076;
1404       break;
1405     case 104792:
1406       fBIFactorA = 0.787215863962;
1407       fBIFactorC = 0.778253173071;
1408       break;
1409     case 104793:
1410       fBIFactorA = 0.692211363661;
1411       fBIFactorC = 0.733152456667;
1412       break;
1413     case 104799:
1414       fBIFactorA = 1.04027825161;
1415       fBIFactorC = 1.00530825942;
1416       break;
1417     case 104800:
1418       fBIFactorA = 1.05309910671;
1419       fBIFactorC = 1.00376801855;
1420       break;
1421     case 104801:
1422       fBIFactorA = 1.0531231922;
1423       fBIFactorC = 0.992439666758;
1424       break;
1425     case 104802:
1426       fBIFactorA = 1.04191478134;
1427       fBIFactorC = 0.979368585208;
1428       break;
1429     case 104803:
1430       fBIFactorA = 1.03121314094;
1431       fBIFactorC = 0.973379962609;
1432       break;
1433     case 104824:
1434       fBIFactorA = 0.969945926722;
1435       fBIFactorC = 0.39549745806;
1436       break;
1437     case 104825:
1438       fBIFactorA = 0.968627213937;
1439       fBIFactorC = 0.310100412205;
1440       break;
1441     case 104841:
1442       fBIFactorA = 0.991601393212;
1443       fBIFactorC = 0.83762204722;
1444       break;
1445     case 104845:
1446       fBIFactorA = 0.98040863886;
1447       fBIFactorC = 0.694824205793;
1448       break;
1449     case 104867:
1450       fBIFactorA = 1.10646173412;
1451       fBIFactorC = 0.841407246916;
1452       break;
1453     case 104876:
1454       fBIFactorA = 1.12063452421;
1455       fBIFactorC = 0.78726542895;
1456       break;
1457     case 104890:
1458       fBIFactorA = 1.02346137453;
1459       fBIFactorC = 1.03355663595;
1460       break;
1461     case 104892:
1462       fBIFactorA = 1.05406025913;
1463       fBIFactorC = 1.00029166135;
1464       break;
1465     case 105143:
1466       fBIFactorA = 0.947343384349;
1467       fBIFactorC = 0.972637444408;
1468       break;
1469     case 105160:
1470       fBIFactorA = 0.908854622177;
1471       fBIFactorC = 0.958851103977;
1472       break;
1473     case 105256:
1474       fBIFactorA = 0.810076150206;
1475       fBIFactorC = 0.884663561883;
1476       break;
1477     case 105257:
1478       fBIFactorA = 0.80974912303;
1479       fBIFactorC = 0.878859123479;
1480       break;
1481     case 105268:
1482       fBIFactorA = 0.809052110679;
1483       fBIFactorC = 0.87233890989;
1484       break;
1485     default:
1486       fBIFactorA = 1;
1487       fBIFactorC = 1;
1488     }
1489   }
1490
1491 }