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