ANALYSIS/AliPhysicsSelection.cxx

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