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