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