ANALYSIS/AliPhysicsSelection.cxx

   1 /* $Id: AliPhysicsSelection.cxx 35782 2009-10-22 11:54:31Z jgrosseo $ */
   2
   3 /**************************************************************************
   4  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   5  *                                                                        *
   6  * Author: The ALICE Off-line Project.                                    *
   7  * Contributors are mentioned in the code where appropriate.              *
   8  *                                                                        *
   9  * Permission to use, copy, modify and distribute this software and its   *
  10  * documentation strictly for non-commercial purposes is hereby granted   *
  11  * without fee, provided that the above copyright notice appears in all   *
  12  * copies and that both the copyright notice and this permission notice   *
  13  * appear in the supporting documentation. The authors make no claims     *
  14  * about the suitability of this software for any purpose. It is          *
  15  * provided "as is" without express or implied warranty.                  *
  16  **************************************************************************/
  17
  18 //-------------------------------------------------------------------------
  19 //                      Implementation of   Class AliPhysicsSelection
  20 // This class selects collision candidates from data runs, applying selection cuts on triggers
  21 // and background rejection based on the content of the ESD
  22 //
  23 // Usage:
  24 //
  25 // Create the object:
  26 //   fPhysicsSelection = new AliPhysicsSelection;
  27 //
  28 // For MC data, call
  29 //   fPhysicsSelection->SetAnalyzeMC()
  30 //
  31 // To check if an event is a collision candidate, use:
  32 //   fPhysicsSelection->IsCollisionCandidate(fESD)
  33 //
  34 // After processing save the resulting histograms to a file with (a folder physics_selection
  35 //   will be created that contains the histograms):
  36 //   fPhysicsSelection->SaveHistograms("physics_selection")
  37 //
  38 // To print statistics after processing use:
  39 //   fPhysicsSelection->Print();
  40 //
  41 // The BX ids corresponding to real bunches crossings p2 are
  42 // automatically selected. You cannot process runs with different
  43 // filling schemes if this option is set. If you want to disable this,
  44 // use:
  45 //   fPhysicsSelection->SetUseBXNumbers(0);
  46 //
  47 //
  48 // If you are analizing muons and you want to keep the muon triggers
  49 // besides the CINT1B you can set:
  50 //   fPhysicsSelection->SetUseMuonTriggers();
  51 //
  52 //
  53 // To compute the Background automatically using the control triggers
  54 // use:
  55 //   fPhysicsSelection->SetComputeBG();
  56 // this will show the value of the Beam Gas, accidentals and good
  57 // events as additional rows in the statistic tables, but it will NOT
  58 // subtract the background automatically.
  59 // This option will only work for runs taken with the CINT1
  60 // suite. This options enables automatically also the usage of BX
  61 // numbers. You can only process one run at a time if you require this
  62 // options, because it uses the bunch intensity estimated run by run.
  63 //
  64 // The BG will usually be more important in the so-called "bin 0": the
  65 // class can also compute the statistics table for events in this
  66 // bin. Since the definition of bin 0 may in general change from
  67 // analysis to analysis, the user needs to provide a callback
  68 // implementing the definition of bin zero. The callback should be
  69 // implemented as a method in the analysis task and should override
  70 // the IsEventInBinZero method of AliAnalysisTaskSE, and should thus
  71 // have the the following prototype:
  72 //   Bool_t IsEventInBinZero();
  73 // It should return true if the event is in the bin 0 and it is set by
  74 // passing to the physics selection the NAME of the task where the
  75 // callback is implemented:
  76 //   fPhysicsSelection->SetBin0Callback("MyTask").
  77 //
  78 //
  79 // Usually the class selects the trigger scheme by itself depending on the run number.
  80 // Nevertheless, you can do that manually by calling AddCollisionTriggerClass() and AddBGTriggerClass()
  81 // Example:
  82 // To define the class CINT1B-ABCE-NOPF-ALL as collision trigger (those will be accepted as
  83 // collision candidates when they pass the selection):
  84 //   AddCollisionTriggerClass("+CINT1B-ABCE-NOPF-ALL #769 #3119");
  85 // To select on bunch crossing IDs in addition, use:
  86 //   AddCollisionTriggerClass("+CINT1B-ABCE-NOPF-ALL #769 #3119");
  87 // To define the class CINT1A-ABCE-NOPF-ALL as a background trigger (those will only be counted
  88 // for the control histograms):
  89 //   AddBGTriggerClass("+CINT1A-ABCE-NOPF-ALL");
  90 // You can also specify more than one trigger class in a string or you can require that some are *not*
  91 // present. The following line would require CSMBA-ABCE-NOPF-ALL, but CSMBB-ABCE-NOPF-ALL is not allowed
  92 // to be present:
  93 //   AddBGTriggerClass("+CSMBA-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL");
  94 //
  95 //   Origin: Jan Fiete Grosse-Oetringhaus, CERN
  96 //           Michele Floris, CERN
  97 //-------------------------------------------------------------------------
  98
  99 #include <Riostream.h>
 100 #include <TH1F.h>
 101 #include <TH2F.h>
 102 #include <TList.h>
 103 #include <TIterator.h>
 104 #include <TDirectory.h>
 105 #include <TObjArray.h>
 106
 107 #include <AliPhysicsSelection.h>
 108
 109 #include <AliTriggerAnalysis.h>
 110 #include <AliLog.h>
 111
 112 #include <AliESDEvent.h>
 113 #include <AliAnalysisTaskSE.h>
 114 #include "AliAnalysisManager.h"
 115 #include "TPRegexp.h"
 116
 117 ClassImp(AliPhysicsSelection)
 118
 119 AliPhysicsSelection::AliPhysicsSelection() :
 120   AliAnalysisCuts("AliPhysicsSelection", "AliPhysicsSelection"),
 121   fCurrentRun(-1),
 122   fMC(kFALSE),
 123   fCollTrigClasses(),
 124   fBGTrigClasses(),
 125   fTriggerAnalysis(),
 126   fBackgroundIdentification(0),
 127   fHistBunchCrossing(0),
 128   fHistTriggerPattern(0),
 129   fSkipTriggerClassSelection(0),
 130   fUsingCustomClasses(0),
 131   fSkipV0(0),
 132   fBIFactorA(1),
 133   fBIFactorC(1),
 134   fBIFactorAC(1),
 135   fComputeBG(0),
 136   fBGStatOffset(0),
 137   fUseBXNumbers(1),
 138   fUseMuonTriggers(0),
 139   fFillingScheme(""),
 140   fBin0CallBack(""),
 141   fBin0CallBackPointer(0)
 142 {
 143   // constructor
 144
 145   fCollTrigClasses.SetOwner(1);
 146   fBGTrigClasses.SetOwner(1);
 147   fTriggerAnalysis.SetOwner(1);
 148   fHistStatistics[0] = 0;
 149   fHistStatistics[1] = 0;
 150
 151   AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kWarning);
 152 }
 153
 154 AliPhysicsSelection::~AliPhysicsSelection()
 155 {
 156   // destructor
 157
 158   fCollTrigClasses.Delete();
 159   fBGTrigClasses.Delete();
 160   fTriggerAnalysis.Delete();
 161
 162   if (fHistStatistics[0])
 163   {
 164     delete fHistStatistics[0];
 165     fHistStatistics[0] = 0;
 166   }
 167   if (fHistStatistics[1])
 168   {
 169     delete fHistStatistics[1];
 170     fHistStatistics[1] = 0;
 171   }
 172
 173   if (fHistBunchCrossing)
 174   {
 175     delete fHistBunchCrossing;
 176     fHistBunchCrossing = 0;
 177   }
 178   if (fHistTriggerPattern)
 179   {
 180     delete fHistTriggerPattern;
 181     fHistTriggerPattern = 0;
 182   }
 183   if (fBackgroundIdentification)
 184   {
 185     delete fBackgroundIdentification;
 186     fBackgroundIdentification = 0;
 187   }
 188 }
 189
 190 UInt_t AliPhysicsSelection::CheckTriggerClass(const AliESDEvent* aEsd, const char* trigger, Int_t& triggerLogic) const
 191 {
 192   // checks if the given trigger class(es) are found for the current event
 193   // format of trigger: +TRIGGER1 -TRIGGER2 [#XXX] [&YY] [*ZZ]
 194   //   requires TRIGGER1 and rejects TRIGGER2
 195   //   in bunch crossing XXX
 196   //   if successful, YY is returned (for association between entry in fCollTrigClasses and AliVEvent::EOfflineTriggerTypes)
 197   //   triggerLogic is filled with ZZ, defaults to kCINT1
 198
 199   Bool_t foundBCRequirement = kFALSE;
 200   Bool_t foundCorrectBC = kFALSE;
 201
 202   UInt_t returnCode = AliVEvent::kUserDefined;
 203   triggerLogic = kCINT1;
 204
 205   TString str(trigger);
 206   TObjArray* tokens = str.Tokenize(" ");
 207
 208   for (Int_t i=0; i < tokens->GetEntries(); i++)
 209   {
 210     TString str2(((TObjString*) tokens->At(i))->String());
 211
 212     if (str2[0] == '+' || str2[0] == '-')
 213     {
 214       Bool_t flag = (str2[0] == '+');
 215
 216       str2.Remove(0, 1);
 217
 218       if (flag && !aEsd->IsTriggerClassFired(str2))
 219       {
 220         AliDebug(AliLog::kDebug, Form("Rejecting event because trigger class %s is not present", str2.Data()));
 221         delete tokens;
 222         return kFALSE;
 223       }
 224       if (!flag && aEsd->IsTriggerClassFired(str2))
 225       {
 226         AliDebug(AliLog::kDebug, Form("Rejecting event because trigger class %s is present", str2.Data()));
 227         delete tokens;
 228         return kFALSE;
 229       }
 230     }
 231     else if (str2[0] == '#')
 232     {
 233       foundBCRequirement = kTRUE;
 234
 235       str2.Remove(0, 1);
 236
 237       Int_t bcNumber = str2.Atoi();
 238       AliDebug(AliLog::kDebug, Form("Checking for bunch crossing number %d", bcNumber));
 239
 240       if (aEsd->GetBunchCrossNumber() == bcNumber)
 241       {
 242         foundCorrectBC = kTRUE;
 243         AliDebug(AliLog::kDebug, Form("Found correct bunch crossing %d", bcNumber));
 244       }
 245     }
 246     else if (str2[0] == '&' && !fUsingCustomClasses)
 247     {
 248       str2.Remove(0, 1);
 249
 250       returnCode = str2.Atoll();
 251     }
 252     else if (str2[0] == '*')
 253     {
 254       str2.Remove(0, 1);
 255
 256       triggerLogic = str2.Atoi();
 257     }
 258     else
 259       AliFatal(Form("Invalid trigger syntax: %s", trigger));
 260   }
 261
 262   delete tokens;
 263
 264   if (foundBCRequirement && !foundCorrectBC)
 265     return kFALSE;
 266
 267   return returnCode;
 268 }
 269
 270 //______________________________________________________________________________
 271 TObject *AliPhysicsSelection::GetStatistics(Option_t *option) const
 272 {
 273 // Get the statistics histograms ("ALL" and "BIN0")
 274    TString opt(option);
 275    opt.ToUpper();
 276    Int_t ihist = 0;
 277    if (opt == "ALL") ihist = kStatIdxAll;
 278    if (opt == "BIN0") ihist = kStatIdxBin0;
 279    return fHistStatistics[ihist];
 280 }
 281
 282 //______________________________________________________________________________
 283 UInt_t AliPhysicsSelection::IsCollisionCandidate(const AliESDEvent* aEsd)
 284 {
 285   // checks if the given event is a collision candidate
 286   //
 287   // returns a bit word describing the fired offline triggers (see AliVEvent::EOfflineTriggerTypes)
 288
 289   if (fCurrentRun != aEsd->GetRunNumber()) {
 290     if (!Initialize(aEsd->GetRunNumber()))
 291       AliFatal(Form("Could not initialize for run %d", aEsd->GetRunNumber()));
 292     if(fComputeBG) SetBIFactors(aEsd); // is this safe here?
 293   }
 294   const AliESDHeader* esdHeader = aEsd->GetHeader();
 295   if (!esdHeader)
 296   {
 297     AliError("ESD Header could not be retrieved");
 298     return kFALSE;
 299   }
 300
 301   // check event type; should be PHYSICS = 7 for data and 0 for MC
 302   if (!fMC)
 303   {
 304     if (esdHeader->GetEventType() != 7)
 305       return kFALSE;
 306   }
 307   else
 308   {
 309     if (esdHeader->GetEventType() != 0)
 310       AliFatal(Form("Invalid event type for MC: %d", esdHeader->GetEventType()));
 311   }
 312
 313   UInt_t accept = 0;
 314
 315   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
 316   for (Int_t i=0; i < count; i++)
 317   {
 318     const char* triggerClass = 0;
 319     if (i < fCollTrigClasses.GetEntries())
 320       triggerClass = ((TObjString*) fCollTrigClasses.At(i))->String();
 321     else
 322       triggerClass = ((TObjString*) fBGTrigClasses.At(i - fCollTrigClasses.GetEntries()))->String();
 323
 324     AliDebug(AliLog::kDebug, Form("Processing trigger class %s", triggerClass));
 325
 326     AliTriggerAnalysis* triggerAnalysis = static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i));
 327
 328     triggerAnalysis->FillTriggerClasses(aEsd);
 329
 330     Int_t triggerLogic = 0;
 331     UInt_t singleTriggerResult = CheckTriggerClass(aEsd, triggerClass, triggerLogic);
 332
 333     if (singleTriggerResult)
 334     {
 335       triggerAnalysis->FillHistograms(aEsd);
 336
 337       Bool_t isBin0 = kFALSE;
 338       if (fBin0CallBack != "") {
 339         AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 340         if (!mgr) {
 341           AliError("Cannot get the analysis manager");
 342         }
 343         else {
 344           isBin0 = ((AliAnalysisTaskSE*)mgr->GetTask(fBin0CallBack.Data()))->IsEventInBinZero();
 345         }
 346       } else if (fBin0CallBackPointer) {
 347           isBin0 = (*fBin0CallBackPointer)(aEsd);
 348       }
 349
 350       // hardware trigger
 351       Int_t fastORHW = triggerAnalysis->SPDFiredChips(aEsd, 1); // SPD number of chips from trigger bits (!)
 352       Int_t fastORHWL1 = triggerAnalysis->SPDFiredChips(aEsd, 1, kFALSE, 2); // SPD number of chips from trigger bits in second layer (!)
 353       Bool_t v0AHW     = fSkipV0 ? 0 :(triggerAnalysis->V0Trigger(aEsd, AliTriggerAnalysis::kASide, kTRUE) == AliTriggerAnalysis::kV0BB);// should replay hw trigger
 354       Bool_t v0CHW     = fSkipV0 ? 0 :(triggerAnalysis->V0Trigger(aEsd, AliTriggerAnalysis::kCSide, kTRUE) == AliTriggerAnalysis::kV0BB);// should replay hw trigger
 355
 356       // offline trigger
 357       Int_t fastOROffline = triggerAnalysis->SPDFiredChips(aEsd, 0); // SPD number of chips from clusters (!)
 358       Int_t fastOROfflineL1 = triggerAnalysis->SPDFiredChips(aEsd, 0, kFALSE, 2); // SPD number of chips from clusters in second layer (!)
 359       Bool_t v0A       = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0A);
 360       Bool_t v0C       = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0C);
 361       Bool_t v0ABG = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0ABG);
 362       Bool_t v0CBG = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0CBG);
 363       Bool_t v0BG = v0ABG || v0CBG;
 364
 365       // fmd
 366       Bool_t fmdA =  triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kFMDA);
 367       Bool_t fmdC =  triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kFMDC);
 368       Bool_t fmd  = fmdA || fmdC;
 369
 370       // SSD
 371       Int_t ssdClusters = triggerAnalysis->SSDClusters(aEsd);
 372
 373       // ZDC
 374       Bool_t zdcA = triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kZDCA);
 375       Bool_t zdcC = triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kZDCC);
 376
 377
 378       // Some "macros"
 379       Bool_t mb1 = (fastOROffline > 0 || v0A || v0C) && (!v0BG);
 380       Bool_t mb1prime = (fastOROffline > 1 || (fastOROffline > 0 && (v0A || v0C)) || (v0A && v0C) ) && (!v0BG);
 381
 382       // Background rejection
 383       Bool_t bgID = kFALSE;
 384       if (fBackgroundIdentification)
 385         bgID = ! fBackgroundIdentification->IsSelected(const_cast<AliESDEvent*> (aEsd));
 386
 387       /*Int_t ntrig = fastOROffline; // any 2 hits
 388       if(v0A)              ntrig += 1;
 389       if(v0C)              ntrig += 1; //v0C alone is enough
 390       if(fmd)              ntrig += 1;
 391       if(ssdClusters>1)    ntrig += 1;*/
 392
 393       // replay hardware trigger (should only remove events for MC)
 394
 395       Bool_t hwTrig = kFALSE;
 396       switch (triggerLogic)
 397       {
 398         case kCINT1:  hwTrig = fastORHW > 0 || v0AHW || v0CHW; break;
 399         case kCMBS2A: hwTrig = fastORHWL1 > 1 && v0AHW; break;
 400         case kCMBS2C: hwTrig = fastORHWL1 > 1 && v0CHW; break;
 401         case kCMBAC:  hwTrig = v0AHW && v0CHW; break;
 402         default: AliFatal(Form("Undefined trigger logic %d", triggerLogic)); break;
 403       }
 404
 405       // Fill trigger pattern histo
 406       Int_t tpatt = 0;
 407       if (fastORHW>0) tpatt+=1;
 408       if (v0AHW)      tpatt+=2;
 409       if (v0CHW)      tpatt+=4;
 410       fHistTriggerPattern->Fill( tpatt );
 411
 412       // fill statistics and return decision
 413       const Int_t nHistStat = 2;
 414       for(Int_t iHistStat = 0; iHistStat < nHistStat; iHistStat++){
 415         if (iHistStat == kStatIdxBin0 && !isBin0) continue; // skip the filling of bin0 stats if the event is not in the bin0
 416
 417         fHistStatistics[iHistStat]->Fill(kStatTriggerClass, i);
 418
 419         // We fill the rest only if hw trigger is ok
 420         if (!hwTrig)
 421           {
 422             AliDebug(AliLog::kDebug, "Rejecting event because hardware trigger is not fired");
 423             continue;
 424           } else {
 425           fHistStatistics[iHistStat]->Fill(kStatHWTrig, i);
 426         }
 427
 428
 429         // v0 BG stats
 430         if (v0ABG)
 431           fHistStatistics[iHistStat]->Fill(kStatV0ABG, i);
 432         if (v0CBG)
 433           fHistStatistics[iHistStat]->Fill(kStatV0CBG, i);
 434
 435         // We fill the rest only if mb1 && ! v0BG
 436         if (mb1)
 437           fHistStatistics[iHistStat]->Fill(kStatMB1, i);
 438         else continue;
 439
 440         if (mb1prime)
 441           fHistStatistics[iHistStat]->Fill(kStatMB1Prime, i);
 442
 443         if (fmd)
 444           fHistStatistics[iHistStat]->Fill(kStatFMD, i);
 445
 446         //if(ntrig >= 2 && !v0BG)
 447         //  fHistStatistics[iHistStat]->Fill(kStatAny2Hits, i);
 448
 449         if (fastOROffline > 0)
 450           fHistStatistics[iHistStat]->Fill(kStatFO1, i);
 451         if (fastOROffline > 1)
 452           fHistStatistics[iHistStat]->Fill(kStatFO2, i);
 453         if (fastOROfflineL1 > 1)
 454           fHistStatistics[iHistStat]->Fill(kStatFO2L1, i);
 455
 456         if (v0A)
 457           fHistStatistics[iHistStat]->Fill(kStatV0A, i);
 458         if (v0C)
 459           fHistStatistics[iHistStat]->Fill(kStatV0C, i);
 460
 461         if (zdcA)
 462           fHistStatistics[iHistStat]->Fill(kStatZDCA, i);
 463         if (zdcC)
 464           fHistStatistics[iHistStat]->Fill(kStatZDCC, i);
 465         if (zdcA && zdcC)
 466           fHistStatistics[iHistStat]->Fill(kStatZDCAC, i);
 467
 468         //       if (fastOROffline > 1 && !v0BG)
 469         //         fHistStatistics[iHistStat]->Fill(kStatFO2NoBG, i);
 470
 471         //if (fastOROffline > 0 && (v0A || v0C) && !v0BG)
 472         //  fHistStatistics[iHistStat]->Fill(kStatFO1AndV0, i);
 473
 474         if (v0A && v0C && !v0BG && !bgID)
 475           fHistStatistics[iHistStat]->Fill(kStatV0, i);
 476
 477         Bool_t offlineAccepted = kFALSE;
 478
 479         switch (triggerLogic)
 480         {
 481           case kCINT1:  offlineAccepted = mb1; break;
 482           case kCMBS2A: offlineAccepted = fastOROfflineL1 > 1 && v0A; break;
 483           case kCMBS2C: offlineAccepted = fastOROfflineL1 > 1 && v0C; break;
 484           case kCMBAC:  offlineAccepted = v0A && v0C; break;
 485           default: AliFatal(Form("Undefined trigger logic %d", triggerLogic)); break;
 486         }
 487
 488         if ( offlineAccepted )
 489           {
 490             if (!v0BG || fSkipV0)
 491               {
 492                 if (!v0BG) fHistStatistics[iHistStat]->Fill(kStatOffline, i);
 493
 494                 if (fBackgroundIdentification && bgID)
 495                   {
 496                     AliDebug(AliLog::kDebug, "Rejecting event because of background identification");
 497                     fHistStatistics[iHistStat]->Fill(kStatBG, i);
 498                   }
 499                 else
 500                   {
 501                     AliDebug(AliLog::kDebug, Form("Accepted event for histograms with trigger logic %d", triggerLogic));
 502
 503                     fHistStatistics[iHistStat]->Fill(kStatAccepted, i);
 504                     if(iHistStat == kStatIdxAll) fHistBunchCrossing->Fill(aEsd->GetBunchCrossNumber(), i); // Fill only for all (avoid double counting)
 505                     if((i < fCollTrigClasses.GetEntries() || fSkipTriggerClassSelection) && (iHistStat==kStatIdxAll))
 506                       accept |= singleTriggerResult; // only set for "all" (should not really matter)
 507                   }
 508               }
 509             else
 510               AliDebug(AliLog::kDebug, "Rejecting event because of V0 BG flag");
 511           }
 512         else
 513           AliDebug(AliLog::kDebug, Form("Rejecting event because trigger logic %d is not fulfilled", triggerLogic));
 514       }
 515     }
 516   }
 517
 518   if (accept)
 519     AliDebug(AliLog::kDebug, Form("Accepted event as collision candidate with bit mask %d", accept));
 520
 521   return accept;
 522 }
 523
 524 Int_t AliPhysicsSelection::GetTriggerScheme(UInt_t runNumber) const
 525 {
 526   // returns the current trigger scheme (classes that are accepted/rejected)
 527
 528   Bool_t pp = kTRUE;
 529   if (runNumber >= 136849)  // last pp run 2010
 530     pp = kFALSE;
 531
 532   if (fMC)
 533   {
 534     if (pp)
 535       return 0; // pp
 536     else
 537       return -1; // HI
 538   }
 539
 540   // TODO dependent on run number
 541
 542   switch (runNumber)
 543   {
 544     // CSMBB triggers
 545     case 104044:
 546     case 105054:
 547     case 105057:
 548       return 2;
 549   }
 550
 551   // defaults
 552   if (pp)
 553     return 1;
 554   else
 555     return 3;
 556 }
 557
 558 const char * AliPhysicsSelection::GetFillingScheme(UInt_t runNumber)  {
 559
 560   if(fMC) return "MC";
 561
 562   if      (runNumber >= 104065 && runNumber <= 104160) {
 563     return "4x4a";
 564   }
 565   else if (runNumber >= 104315 && runNumber <= 104321) {
 566     return "4x4a*";
 567   }
 568   else if (runNumber >= 104792 && runNumber <= 104803) {
 569     return "4x4b";
 570   }
 571   else if (runNumber >= 104824 && runNumber <= 104892) {
 572     return "4x4c";
 573   }
 574   else if (runNumber == 105143 || runNumber == 105160) {
 575     return "16x16a";
 576   }
 577   else if (runNumber >= 105256 && runNumber <= 105268) {
 578     return "4x4c";
 579   }
 580   else if (runNumber >= 114786 && runNumber <= 116684) {
 581     return "Single_2b_1_1_1";
 582   }
 583   else if (runNumber >= 117048 && runNumber <= 117120) {
 584     return "Single_3b_2_2_2";
 585   }
 586   else if (runNumber >= 117220 && runNumber <= 119163) {
 587     return "Single_2b_1_1_1";
 588   }
 589   else if (runNumber >= 119837 && runNumber <= 119862) {
 590     return "Single_4b_2_2_2";
 591   }
 592   else if (runNumber >= 119902 && runNumber <= 120691) {
 593     return "Single_6b_3_3_3";
 594   }
 595   else if (runNumber >= 120741 && runNumber <= 122375) {
 596     return "Single_13b_8_8_8";
 597   }
 598   else if (runNumber >= 130148 && runNumber <= 130375) {
 599     return "125n_48b_36_16_36";
 600   }
 601   else if (runNumber >= 130601 && runNumber <= 130640) {
 602     return "1000ns_50b_35_14_35";
 603   }
 604   else {
 605     AliError(Form("Unknown filling scheme (run %d)", runNumber));
 606   }
 607
 608   return "Unknown";
 609 }
 610
 611 const char * AliPhysicsSelection::GetBXIDs(UInt_t runNumber, const char * trigger)  {
 612
 613   if (!fUseBXNumbers || fMC) return "";
 614
 615   if      (runNumber >= 104065 && runNumber <= 104160) {
 616     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #2128 #3019";
 617     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #346 #3465";
 618     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1234 #1680";
 619     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
 620     //    else AliError(Form("Unknown trigger: %s", trigger));
 621   }
 622   else if (runNumber >= 104315 && runNumber <= 104321) {
 623     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #2000 #2891";
 624     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #218 #3337";
 625     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1106 #1552";
 626     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
 627     //    else AliError(Form("Unknown trigger: %s", trigger));
 628   }
 629   else if (runNumber >= 104792 && runNumber <= 104803) {
 630     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #2228 #3119";
 631     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2554 #446";
 632     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1334 #769";
 633     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
 634     //    else AliError(Form("Unknown trigger: %s", trigger));
 635   }
 636   else if (runNumber >= 104824 && runNumber <= 104892) {
 637     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #3119 #769";
 638     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2554 #446";
 639     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1334 #2228";
 640     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
 641     //    else AliError(Form("Unknown trigger: %s", trigger));
 642   }
 643   else if (runNumber == 105143 || runNumber == 105160) {
 644     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #1337 #1418 #2228 #2309 #3119 #3200 #446 #527";
 645     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #1580  #1742  #1904  #2066  #2630  #2792  #2954  #3362";
 646     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "  #845  #1007  #1169   #1577 #3359 #3521 #119  #281 ";
 647     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
 648     //    else AliError(Form("Unknown trigger: %s", trigger));
 649   }
 650   else if (runNumber >= 105256 && runNumber <= 105268) {
 651     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #3019 #669";
 652     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2454 #346";
 653     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1234 #2128";
 654     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
 655     //    else AliError(Form("Unknown trigger: %s", trigger));
 656   } else if (runNumber >= 114786 && runNumber <= 116684) { // 7 TeV 2010, assume always the same filling scheme
 657     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #346";
 658     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2131";
 659     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #3019";
 660     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #1238";
 661     //    else AliError(Form("Unknown trigger: %s", trigger));
 662   }
 663   else if (runNumber >= 117048 && runNumber <= 117120) {
 664     //    return "Single_3b_2_2_2";
 665    if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346  #1240 ";
 666    else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #2131 ";
 667    else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019 ";
 668    else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1238";
 669    //   else AliError(Form("Unknown trigger: %s", trigger));
 670
 671   }
 672   else if ((runNumber >= 117220 && runNumber <= 118555) || (runNumber >= 118784 && runNumber <= 119163))
 673   {
 674     //    return "Single_2b_1_1_1";
 675     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346 ";
 676     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #2131 ";
 677     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019 ";
 678     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1238 ";
 679     //    else AliError(Form("Unknown trigger: %s", trigger));
 680   }
 681   else if (runNumber >= 118556 && runNumber <= 118783) {
 682     //    return "Single_2b_1_1_1";
 683     // same as previous but was misaligned by 1 BX in fill 1069
 684     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #345 ";
 685     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #2130 ";
 686     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3018 ";
 687     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1238 ";
 688     //    else AliError(Form("Unknown trigger: %s", trigger));
 689   }
 690   else if (runNumber >= 119837 && runNumber <= 119862) {
 691     //    return "Single_4b_2_2_2";
 692     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #669  #3019 ";
 693     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #346  #2454 ";
 694     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #1234  #2128 ";
 695     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1681 #3463";
 696     //    else AliError(Form("Unknown trigger: %s", trigger));
 697
 698   }
 699   else if (runNumber >= 119902 && runNumber <= 120691) {
 700     //    return "Single_6b_3_3_3";
 701     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346  #546  #746 ";
 702     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #2131  #2331  #2531 ";
 703     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019  #3219  #3419 ";
 704     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1296 #1670";
 705     //    else AliError(Form("Unknown trigger: %s", trigger));
 706   }
 707   else if (runNumber >= 120741 && runNumber <= 122375) {
 708     //    return "Single_13b_8_8_8";
 709     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346  #446  #546  #646  #1240  #1340  #1440  #1540 ";
 710     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #946  #2131  #2231  #2331  #2431 ";
 711     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019  #3119  #3219  #3319  #3519 ";
 712     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1835 #2726";
 713     //    else AliError(Form("Unknown trigger: %s", trigger));
 714
 715   }
 716   else if (runNumber >= 130148 && runNumber <= 130375) {
 717     TString triggerString = trigger;
 718     static TString returnString = " ";
 719     returnString = "";
 720     if (triggerString.Contains("B")) returnString += "   #346  #396  #446  #496  #546  #596  #646  #696  #1240  #1290  #1340  #1390  #1440  #1490  #1540  #1590 ";
 721     if (triggerString.Contains("A")) returnString += "   #755  #805  #855  #905  #955  #1005  #1799  #1849  #1899  #2131  #2181  #2231  #2281  #2331  #2381  #2431  #2481  #2531  #2581  #2631  #2846  #3016  #3066  #3116  #3166  #3216  #3266  #3316  #3366  #3425  #3475  #3525 ";
 722     if (triggerString.Contains("C")) returnString += "   #3019  #3069  #3119  #3169  #3219  #3269  #3319  #3369  #14  #64  #114  #746  #796  #846  #908  #958  #1008  #1640  #1690  #1740  #2055  #2125  #2175  #2225  #2275  #2325  #2375  #2425  #2475  #2534  #2584  #2634 ";
 723     // Printf("0x%x",returnString.Data());
 724     // Printf("%s",returnString.Data());
 725     return returnString.Data();
 726   }
 727   else if (runNumber >= 130601 && runNumber <= 130640) {
 728     TString triggerString = trigger;
 729     static TString returnString = " ";
 730     returnString = "";
 731     if (triggerString.Contains("B")) returnString += "  #346  #386  #426  #466  #506  #546  #586  #1240  #1280  #1320  #1360  #1400  #1440  #1480 ";
 732     if (triggerString.Contains("A")) returnString += "  #626  #666  #706  #746  #786  #826  #866  #1520  #1560  #1600  #1640  #1680  #1720  #1760  #2076  #2131  #2171  #2211  #2251  #2291  #2331  #2371  #2414  #2454  #2494  #2534  #2574  #2614  #2654  #2694  #2734  #2774  #2814 "; //#2854  #2894  #2934 not present in this run
 733     if (triggerString.Contains("C")) returnString += "  #3019  #3059  #3099  #3139  #3179  #3219  #3259  #3299  #3339  #3379  #3419  #3459  #3499  #3539  #115  #629  #669  #709  #749  #789  #829  #869  #909  #949  #989  #1029  #1069  #1109  #1149  #1523  #1563  #1603  #1643 "; //#1683  #1723  #1763 not present in this run
 734     return returnString.Data();
 735   }
 736
 737   else {
 738     AliWarning(Form("Unknown run %d, using all BXs!",runNumber));
 739   }
 740
 741   return "";
 742 }
 743
 744 Bool_t AliPhysicsSelection::Initialize(Int_t runNumber)
 745 {
 746   // initializes the object for the given run
 747   // 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
 748
 749   Bool_t oldStatus = TH1::AddDirectoryStatus();
 750   TH1::AddDirectory(kFALSE);
 751
 752   if(!fBin0CallBack)
 753     AliError("Bin0 Callback not set: will not fill the statistics for the bin 0");
 754
 755   if (fMC) {
 756     // override BX and bg options in case of MC
 757     fComputeBG    = kFALSE;
 758     fUseBXNumbers = kFALSE;
 759   }
 760
 761   Int_t triggerScheme = GetTriggerScheme(runNumber);
 762   if (!fUsingCustomClasses && fCurrentRun != -1 && triggerScheme != GetTriggerScheme(fCurrentRun))
 763     AliFatal("Processing several runs with different trigger schemes is not supported");
 764
 765   if(fComputeBG && fCurrentRun != -1 && fCurrentRun != runNumber)
 766     AliFatal("Cannot process several runs because BG computation is requested");
 767
 768   if(fComputeBG && !fUseBXNumbers)
 769     AliFatal("Cannot compute BG if BX numbers are not used");
 770
 771   if(fUseBXNumbers && fFillingScheme != "" && fFillingScheme != GetFillingScheme(runNumber))
 772     AliFatal("Cannot process runs with different filling scheme if usage of BX numbers is requested");
 773
 774   fFillingScheme      = GetFillingScheme(runNumber);
 775
 776   AliInfo(Form("Initializing for run %d", runNumber));
 777
 778   // initialize first time?
 779   if (fCurrentRun == -1)
 780   {
 781     if (fUsingCustomClasses) {
 782       AliInfo("Using user-provided trigger classes");
 783     } else {
 784       switch (triggerScheme)
 785       {
 786       case -1:
 787         // MC Heavy-Ion
 788
 789         fCollTrigClasses.Add(new TObjString(Form("&%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2A)));
 790         fCollTrigClasses.Add(new TObjString(Form("&%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2C)));
 791         fCollTrigClasses.Add(new TObjString(Form("&%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBAC)));
 792         break;
 793
 794       case 0:
 795         // MC Proton-Proton
 796
 797         fCollTrigClasses.Add(new TObjString(Form("&%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCINT1)));
 798         break;
 799
 800       case 1:
 801         // Proton-Proton
 802
 803         // trigger classes used before August 2010
 804         fCollTrigClasses.Add(new TObjString(Form("%s%s &%u","+CINT1B-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1B-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
 805         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CINT1A-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1A-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
 806         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CINT1C-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1C-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
 807         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CINT1-E-NOPF-ALL",      GetBXIDs(runNumber,"CINT1-E-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
 808
 809         // Muon trigger have the same BXIDs of the corresponding CINT triggers
 810         fCollTrigClasses.Add(new TObjString(Form("%s%s &%u","+CMUS1B-ABCE-NOPF-MUON",  GetBXIDs(runNumber,"CINT1B-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
 811         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CMUS1A-ABCE-NOPF-MUON",  GetBXIDs(runNumber,"CINT1A-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
 812         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CMUS1C-ABCE-NOPF-MUON",  GetBXIDs(runNumber,"CINT1C-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
 813         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CMUS1-E-NOPF-MUON"    ,  GetBXIDs(runNumber,"CINT1-E-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
 814
 815         // triggers classes used from August 2010
 816         // MB
 817         fCollTrigClasses.Add(new TObjString(Form("%s%s &%u", "+CINT1-B-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "B"),  (UInt_t) AliVEvent::kMB)));
 818         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CINT1-AC-NOPF-ALLNOTRD", GetBXIDs(runNumber, "AC"), (UInt_t) AliVEvent::kMB)));
 819         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CINT1-E-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "E"),  (UInt_t) AliVEvent::kMB)));
 820
 821         // MUON
 822         fCollTrigClasses.Add(new TObjString(Form("%s%s &%u", "+CMUS1-B-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "B"),  (UInt_t) AliVEvent::kMUON)));
 823         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CMUS1-AC-NOPF-ALLNOTRD", GetBXIDs(runNumber, "AC"), (UInt_t) AliVEvent::kMUON)));
 824         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CMUS1-E-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "E"),  (UInt_t) AliVEvent::kMUON)));
 825
 826         // High Multiplicity
 827         fCollTrigClasses.Add(new TObjString(Form("%s%s &%u", "+CSH1-B-NOPF-ALLNOTRD"  , GetBXIDs(runNumber, "B"),  (UInt_t) AliVEvent::kHighMult)));
 828         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CSH1-AC-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "AC"), (UInt_t) AliVEvent::kHighMult)));
 829         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CSH1-E-NOPF-ALLNOTRD"  , GetBXIDs(runNumber, "E"),  (UInt_t) AliVEvent::kHighMult)));
 830
 831         // WARNING: IF YOU ADD MORE TRIGGER CLASSES, PLEASE CHECK THAT THE REGULAR EXPRESSION IN GetStatRow IS STILL VALID
 832
 833         break;
 834
 835       case 2:
 836         // Proton-Proton
 837
 838         fCollTrigClasses.Add(new TObjString(Form("+CSMBB-ABCE-NOPF-ALL &%u", (UInt_t) AliVEvent::kMB)));
 839         fBGTrigClasses.Add(new TObjString(Form("+CSMBA-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL &%u", (UInt_t) AliVEvent::kMB)));
 840         fBGTrigClasses.Add(new TObjString(Form("+CSMBC-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL &%u", (UInt_t) AliVEvent::kMB)));
 841         break;
 842
 843       case 3:
 844         // Heavy-ion
 845
 846         fCollTrigClasses.Add(new TObjString(Form("+CMBS2A-B-NOPF-ALL &%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2A)));
 847         fBGTrigClasses.Add  (new TObjString(Form("+CMBS2A-A-NOPF-ALL &%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2A)));
 848         fBGTrigClasses.Add  (new TObjString(Form("+CMBS2A-C-NOPF-ALL &%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2A)));
 849         fBGTrigClasses.Add  (new TObjString(Form("+CMBS2A-E-NOPF-ALL &%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2A)));
 850
 851         fCollTrigClasses.Add(new TObjString(Form("+CMBS2C-B-NOPF-ALL &%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2C)));
 852         fBGTrigClasses.Add  (new TObjString(Form("+CMBS2C-A-NOPF-ALL &%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2C)));
 853         fBGTrigClasses.Add  (new TObjString(Form("+CMBS2C-C-NOPF-ALL &%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2C)));
 854         fBGTrigClasses.Add  (new TObjString(Form("+CMBS2C-E-NOPF-ALL &%u *%d", (UInt_t) AliVEvent::kMB, (Int_t) kCMBS2C)));
 855
 856         fCollTrigClasses.Add(new TObjString(Form("+CMBAC-B-NOPF-ALL &%u *%d",  (UInt_t) AliVEvent::kMB, (Int_t) kCMBAC)));
 857         fBGTrigClasses.Add  (new TObjString(Form("+CMBAC-A-NOPF-ALL &%u *%d",  (UInt_t) AliVEvent::kMB, (Int_t) kCMBAC)));
 858         fBGTrigClasses.Add  (new TObjString(Form("+CMBAC-C-NOPF-ALL &%u *%d",  (UInt_t) AliVEvent::kMB, (Int_t) kCMBAC)));
 859         fBGTrigClasses.Add  (new TObjString(Form("+CMBAC-E-NOPF-ALL &%u *%d",  (UInt_t) AliVEvent::kMB, (Int_t) kCMBAC)));
 860
 861         break;
 862
 863       default:
 864         AliFatal(Form("Unsupported trigger scheme %d", triggerScheme));
 865       }
 866     }
 867
 868     Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
 869
 870     for (Int_t i=0; i<count; i++)
 871     {
 872       AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis;
 873       triggerAnalysis->SetAnalyzeMC(fMC);
 874       triggerAnalysis->EnableHistograms();
 875       triggerAnalysis->SetSPDGFOThreshhold(1);
 876       triggerAnalysis->SetDoFMD(kFALSE);
 877       fTriggerAnalysis.Add(triggerAnalysis);
 878     }
 879
 880     // TODO: shall I really delete this?
 881     if (fHistStatistics[0])
 882       delete fHistStatistics[0];
 883     if (fHistStatistics[1])
 884       delete fHistStatistics[1];
 885
 886     fHistStatistics[kStatIdxBin0] = BookHistStatistics("_Bin0");
 887     fHistStatistics[kStatIdxAll]  = BookHistStatistics("");
 888
 889     if (fHistBunchCrossing)
 890       delete fHistBunchCrossing;
 891
 892     fHistBunchCrossing = new TH2F("fHistBunchCrossing", "fHistBunchCrossing;bunch crossing number;", 4000, -0.5, 3999.5,  count, -0.5, -0.5 + count);
 893
 894     if (fHistTriggerPattern)
 895       delete fHistTriggerPattern;
 896
 897     const int ntrig=3;
 898     Int_t n = 1;
 899     const Int_t nbinTrig = TMath::Nint(TMath::Power(2,ntrig));
 900
 901     fHistTriggerPattern = new TH1F("fHistTriggerPattern", "Trigger pattern: FO + 2*v0A + 4*v0C",
 902                                    nbinTrig, -0.5, nbinTrig-0.5);
 903     fHistTriggerPattern->GetXaxis()->SetBinLabel(1,"NO TRIG");
 904     fHistTriggerPattern->GetXaxis()->SetBinLabel(2,"FO");
 905     fHistTriggerPattern->GetXaxis()->SetBinLabel(3,"v0A");
 906     fHistTriggerPattern->GetXaxis()->SetBinLabel(4,"FO & v0A");
 907     fHistTriggerPattern->GetXaxis()->SetBinLabel(5,"v0C");
 908     fHistTriggerPattern->GetXaxis()->SetBinLabel(6,"FO & v0C");
 909     fHistTriggerPattern->GetXaxis()->SetBinLabel(7,"v0A & v0C");
 910     fHistTriggerPattern->GetXaxis()->SetBinLabel(8,"FO & v0A & v0C");
 911
 912
 913     n = 1;
 914     for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
 915     {
 916       fHistBunchCrossing->GetYaxis()->SetBinLabel(n, ((TObjString*) fCollTrigClasses.At(i))->String());
 917       n++;
 918     }
 919     for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
 920     {
 921       fHistBunchCrossing->GetYaxis()->SetBinLabel(n, ((TObjString*) fBGTrigClasses.At(i))->String());
 922       n++;
 923     }
 924
 925
 926
 927   }
 928
 929   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
 930   for (Int_t i=0; i<count; i++)
 931   {
 932     AliTriggerAnalysis* triggerAnalysis = static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i));
 933
 934     switch (runNumber)
 935     {
 936       case 104315:
 937       case 104316:
 938       case 104320:
 939       case 104321:
 940       case 104439:
 941         triggerAnalysis->SetV0TimeOffset(7.5);
 942         break;
 943       default:
 944         triggerAnalysis->SetV0TimeOffset(0);
 945     }
 946   }
 947
 948   fCurrentRun = runNumber;
 949
 950   TH1::AddDirectory(oldStatus);
 951
 952   return kTRUE;
 953 }
 954
 955 TH2F * AliPhysicsSelection::BookHistStatistics(const char * tag) {
 956   // add 6 rows to count for the estimate of good, accidentals and
 957   // BG and the ratio of BG and accidentals to total +ratio goot to
 958   // first col + 2 for error on good.
 959   // TODO: Remember the the indexes of rows for the BG selection. Add new member fBGRows[] and use kStat as indexes
 960
 961   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
 962 #ifdef VERBOSE_STAT
 963   Int_t extrarows = fComputeBG != 0 ? 8 : 0;
 964 #else
 965   Int_t extrarows = fComputeBG != 0 ? 3 : 0;
 966 #endif
 967   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);
 968
 969   h->GetXaxis()->SetBinLabel(kStatTriggerClass,  "Trigger class");
 970   h->GetXaxis()->SetBinLabel(kStatHWTrig,        "Hardware trigger");
 971   h->GetXaxis()->SetBinLabel(kStatFO1,           "FO >= 1");
 972   h->GetXaxis()->SetBinLabel(kStatFO2,           "FO >= 2");
 973   h->GetXaxis()->SetBinLabel(kStatFO2L1,         "FO (L1) >= 2");
 974   h->GetXaxis()->SetBinLabel(kStatV0A,           "V0A");
 975   h->GetXaxis()->SetBinLabel(kStatV0C,           "V0C");
 976   h->GetXaxis()->SetBinLabel(kStatFMD,           "FMD");
 977   h->GetXaxis()->SetBinLabel(kStatV0ABG,         "V0A BG");
 978   h->GetXaxis()->SetBinLabel(kStatV0CBG,         "V0C BG");
 979   h->GetXaxis()->SetBinLabel(kStatZDCA,          "ZDCA");
 980   h->GetXaxis()->SetBinLabel(kStatZDCC,          "ZDCC");
 981   h->GetXaxis()->SetBinLabel(kStatZDCAC,         "ZDCA & ZDCC");
 982   h->GetXaxis()->SetBinLabel(kStatMB1,           "(FO >= 1 | V0A | V0C) & !V0 BG");
 983   h->GetXaxis()->SetBinLabel(kStatMB1Prime,      "(FO >= 2 | (FO >= 1 & (V0A | V0C)) | (V0A &v0C) ) & !V0 BG");
 984   //h->GetXaxis()->SetBinLabel(kStatFO1AndV0,    "FO >= 1 & (V0A | V0C) & !V0 BG");
 985   h->GetXaxis()->SetBinLabel(kStatV0,            "V0A & V0C & !V0 BG & !BG ID");
 986   h->GetXaxis()->SetBinLabel(kStatOffline,       "Offline Trigger");
 987   //h->GetXaxis()->SetBinLabel(kStatAny2Hits,    "2 Hits & !V0 BG");
 988   h->GetXaxis()->SetBinLabel(kStatBG,            "Background identification");
 989   h->GetXaxis()->SetBinLabel(kStatAccepted,      "Accepted");
 990
 991   Int_t n = 1;
 992   for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
 993     {
 994       h->GetYaxis()->SetBinLabel(n, ((TObjString*) fCollTrigClasses.At(i))->String());
 995       n++;
 996     }
 997   for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
 998     {
 999       h->GetYaxis()->SetBinLabel(n, ((TObjString*) fBGTrigClasses.At(i))->String());
1000       n++;
1001     }
1002
1003   if(fComputeBG) {
1004     fBGStatOffset = n;
1005     h->GetYaxis()->SetBinLabel(n++, Form("BG (A+C) (Mask [0x%x])", fComputeBG));
1006     h->GetYaxis()->SetBinLabel(n++, "ACC");
1007 #ifdef VERBOSE_STAT
1008     h->GetYaxis()->SetBinLabel(n++, "BG (A+C) %  (rel. to CINT1B)");
1009     h->GetYaxis()->SetBinLabel(n++, "ACC % (rel. to CINT1B)");
1010     h->GetYaxis()->SetBinLabel(n++, "ERR GOOD %");
1011     h->GetYaxis()->SetBinLabel(n++, "GOOD % (rel. to 1st col)");
1012     h->GetYaxis()->SetBinLabel(n++, "ERR GOOD");
1013 #endif
1014     h->GetYaxis()->SetBinLabel(n++, "GOOD");
1015   }
1016
1017   return h;
1018 }
1019
1020 void AliPhysicsSelection::Print(Option_t *option) const
1021 {
1022   // print the configuration
1023   TString msg;
1024   Printf("Configuration initialized for run %d (MC: %d):", fCurrentRun, fMC);
1025   msg += Form("Configuration initialized for run %d (MC: %d):\n", fCurrentRun, fMC);
1026
1027   Printf("Collision trigger classes:");
1028   for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
1029     Printf("%s", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
1030
1031   Printf("Background trigger classes:");
1032   for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
1033     Printf("%s", ((TObjString*) fBGTrigClasses.At(i))->String().Data());
1034
1035   AliTriggerAnalysis* triggerAnalysis = dynamic_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(0));
1036
1037   if (triggerAnalysis)
1038   {
1039     if (triggerAnalysis->GetV0TimeOffset() > 0)
1040       Printf("V0 time offset active: %.2f ns", triggerAnalysis->GetV0TimeOffset());
1041
1042     Printf("\nTotal available events:");
1043
1044     triggerAnalysis->PrintTriggerClasses();
1045   }
1046
1047   if (fHistStatistics[kStatIdxAll])
1048   {
1049     for (Int_t i=0; i<fCollTrigClasses.GetEntries(); i++)
1050     {
1051       Printf("\nSelection statistics for collision trigger %s:", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
1052       msg += Form("\nSelection statistics for collision trigger %s:\n", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
1053
1054       Printf("Total events with correct trigger class: %d", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(1, i+1));
1055       msg += Form("Total events with correct trigger class: %d\n", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(1, i+1));
1056       Printf("Selected collision candidates: %d", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(fHistStatistics[kStatIdxAll]->GetXaxis()->FindBin("Accepted"), i+1));
1057       msg += Form("Selected collision candidates: %d\n", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(fHistStatistics[kStatIdxAll]->GetXaxis()->FindBin("Accepted"), i+1));
1058     }
1059   }
1060
1061   if (fHistBunchCrossing)
1062   {
1063     Printf("\nBunch crossing statistics:");
1064     msg += "\nBunch crossing statistics:\n";
1065
1066     for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
1067     {
1068       TString str;
1069       str.Form("Trigger %s has accepted events in the bunch crossings: ", fHistBunchCrossing->GetYaxis()->GetBinLabel(i));
1070
1071       for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
1072         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
1073           str += Form("%d, ", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
1074
1075       Printf("%s", str.Data());
1076       msg += str;
1077       msg += "\n";
1078     }
1079
1080     for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
1081     {
1082       Int_t countColl = 0;
1083       Int_t countBG = 0;
1084       for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
1085       {
1086         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
1087         {
1088           if (fCollTrigClasses.FindObject(fHistBunchCrossing->GetYaxis()->GetBinLabel(i)))
1089             countColl++;
1090           if (fBGTrigClasses.FindObject(fHistBunchCrossing->GetYaxis()->GetBinLabel(i)))
1091             countBG++;
1092         }
1093       }
1094       if (countColl > 0 && countBG > 0)
1095         Printf("WARNING: Bunch crossing %d has collision and BG trigger classes active. Check BPTX functioning for this run!", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
1096     }
1097   }
1098
1099   if (fUsingCustomClasses)
1100     Printf("WARNING: Using custom trigger classes!");
1101   if (fSkipTriggerClassSelection)
1102     Printf("WARNING: Skipping trigger class selection!");
1103   if (fSkipV0)
1104     Printf("WARNING: Ignoring V0 information in selection");
1105   if(!fBin0CallBack)
1106     Printf("WARNING: Callback not set: will not fill the statistics for the bin 0");
1107   TString opt(option);
1108   opt.ToUpper();
1109   if (opt == "STAT") {
1110      AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1111      if (mgr) mgr->AddStatisticsMsg(msg);
1112   }
1113 }
1114
1115 Long64_t AliPhysicsSelection::Merge(TCollection* list)
1116 {
1117   // Merge a list of AliMultiplicityCorrection objects with this (needed for
1118   // PROOF).
1119   // Returns the number of merged objects (including this).
1120
1121   if (!list)
1122     return 0;
1123
1124   if (list->IsEmpty())
1125     return 1;
1126
1127   TIterator* iter = list->MakeIterator();
1128   TObject* obj;
1129
1130   // collections of all histograms
1131   const Int_t nHists = 9;
1132   TList collections[nHists];
1133
1134   Int_t count = 0;
1135   while ((obj = iter->Next())) {
1136
1137     AliPhysicsSelection* entry = dynamic_cast<AliPhysicsSelection*> (obj);
1138     if (entry == 0)
1139       continue;
1140     // Update run number. If this one is not initialized (-1) take the one from
1141     // the next physics selection to be merged with. In case of 2 different run
1142     // numbers issue a warning (should physics selections from different runs be
1143     // merged together) A.G.
1144     Int_t currentRun = entry->GetCurrentRun();
1145     // Nothing to merge with since run number was not initialized.
1146     if (currentRun < 0) continue;
1147     if (fCurrentRun < 0)
1148     {
1149       fCurrentRun = currentRun;
1150       fMC = entry->fMC;
1151       for (Int_t i=0; i<entry->fCollTrigClasses.GetEntries(); i++)
1152         fCollTrigClasses.Add(entry->fCollTrigClasses.At(i)->Clone());
1153       for (Int_t i=0; i<entry->fBGTrigClasses.GetEntries(); i++)
1154         fBGTrigClasses.Add(entry->fBGTrigClasses.At(i)->Clone());
1155     }
1156     if (fCurrentRun != currentRun)
1157        AliWarning(Form("Current run %d not matching the one to be merged with %d", fCurrentRun, currentRun));
1158
1159     if (entry->fTriggerAnalysis.GetEntries() > 0)
1160       collections[0].Add(&(entry->fTriggerAnalysis));
1161     if (entry->fHistStatistics[0])
1162       collections[1].Add(entry->fHistStatistics[0]);
1163     if (entry->fHistStatistics[1])
1164       collections[2].Add(entry->fHistStatistics[1]);
1165     if (entry->fHistBunchCrossing)
1166       collections[3].Add(entry->fHistBunchCrossing);
1167     if (entry->fHistTriggerPattern)
1168       collections[4].Add(entry->fHistTriggerPattern);
1169     if (entry->fBackgroundIdentification)
1170       collections[5].Add(entry->fBackgroundIdentification);
1171
1172     count++;
1173   }
1174
1175   if (fTriggerAnalysis.GetEntries() == 0 && collections[0].GetEntries() > 0)
1176   {
1177     TList* firstList = (TList*) collections[0].First();
1178     for (Int_t i=0; i<firstList->GetEntries(); i++)
1179       fTriggerAnalysis.Add(firstList->At(i)->Clone());
1180
1181     collections[0].RemoveAt(0);
1182   }
1183   fTriggerAnalysis.Merge(&collections[0]);
1184
1185   // if this instance is empty (not initialized) nothing would be merged here --> copy first entry
1186   if (!fHistStatistics[0] && collections[1].GetEntries() > 0)
1187   {
1188     fHistStatistics[0] = (TH2F*) collections[1].First()->Clone();
1189     collections[1].RemoveAt(0);
1190   }
1191   if (fHistStatistics[0])
1192     fHistStatistics[0]->Merge(&collections[1]);
1193
1194   if (!fHistStatistics[1] && collections[2].GetEntries() > 0)
1195   {
1196     fHistStatistics[1] = (TH2F*) collections[2].First()->Clone();
1197     collections[2].RemoveAt(0);
1198   }
1199   if (fHistStatistics[1])
1200     fHistStatistics[1]->Merge(&collections[2]);
1201
1202   if (!fHistBunchCrossing && collections[3].GetEntries() > 0)
1203   {
1204     fHistBunchCrossing = (TH2F*) collections[3].First()->Clone();
1205     collections[3].RemoveAt(0);
1206   }
1207   if (fHistBunchCrossing)
1208     fHistBunchCrossing->Merge(&collections[3]);
1209
1210   if (!fHistTriggerPattern && collections[4].GetEntries() > 0)
1211   {
1212     fHistTriggerPattern = (TH1F*) collections[4].First()->Clone();
1213     collections[4].RemoveAt(0);
1214   }
1215   if (fHistTriggerPattern)
1216     fHistTriggerPattern->Merge(&collections[4]);
1217
1218   if (!fBackgroundIdentification && collections[5].GetEntries() > 0)
1219   {
1220     fBackgroundIdentification = (AliAnalysisCuts*) collections[5].First()->Clone();
1221     collections[5].RemoveAt(0);
1222   }
1223   if (fBackgroundIdentification)
1224     fBackgroundIdentification->Merge(&collections[5]);
1225
1226   delete iter;
1227
1228   return count+1;
1229 }
1230
1231 void AliPhysicsSelection::SaveHistograms(const char* folder) const
1232 {
1233   // write histograms to current directory
1234
1235   if (!fHistStatistics[0] || !fHistStatistics[1])
1236     return;
1237
1238   if (folder)
1239   {
1240     gDirectory->mkdir(folder);
1241     gDirectory->cd(folder);
1242   }
1243
1244
1245   // Fill the last rows of fHistStatistics before saving
1246   if (fComputeBG) {
1247     Int_t triggerScheme = GetTriggerScheme(UInt_t(fCurrentRun));
1248     if(triggerScheme != 1){
1249       AliWarning("BG estimate only supported for trigger scheme \"1\" (CINT1 suite)");
1250     } else if (fUseMuonTriggers) {
1251       AliWarning("BG estimate with muon triggers to be implemented");
1252     } else {
1253       AliInfo("BG estimate assumes that for a given run you only have A and C triggers separately or"
1254               " toghether as a AC class! Make sure this assumption holds in your case");
1255
1256       // use an anum for the different trigger classes, to make loops easier to read
1257       enum {kClassB =0 , kClassA, kClassC, kClassAC, kClassE, kNClasses};
1258       const char * classFlags[] = {"B", "A", "C", "AC", "E"}; // labels
1259
1260       UInt_t * rows[kNClasses] = {0}; // Array of matching rows
1261       Int_t nrows[kNClasses] = {0};
1262       // Get rows matching the requested trigger bits for all trigger classes
1263       for(Int_t iTrigClass = 0; iTrigClass < kNClasses; iTrigClass++){
1264         nrows[iTrigClass] = GetStatRow(classFlags[iTrigClass],fComputeBG,&rows[iTrigClass]);
1265       }
1266
1267       // 0. Determine the ratios of triggers E/B, A/B, C/B from the stat histogram
1268       // Those are used to rescale the different classes to the same number of bx ids
1269       // TODO: pass names of the rows for B, CA and E and look names of the rows. How do I handle the case in which both AC are in the same row?
1270       Int_t nBXIds[kNClasses] = {0};
1271       cout <<"Computing BG:" << endl;
1272
1273       for(Int_t iTrigClass = 0; iTrigClass < kNClasses; iTrigClass++){
1274         for(Int_t irow = 0; irow < nrows[iTrigClass]; irow++) {
1275           if(irow==0) cout << "- Class " << classFlags[iTrigClass] << endl;
1276           for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++) {
1277             if (fHistBunchCrossing->GetBinContent(j, rows[iTrigClass][irow]) > 0) {
1278               nBXIds[iTrigClass]++;
1279             }
1280           }
1281           if(nBXIds[iTrigClass]>0) cout << "   Using row " << rows[iTrigClass][irow] <<  ": "
1282                                         << fHistBunchCrossing->GetYaxis()->GetBinLabel(rows[iTrigClass][irow])
1283                                         << " (nBXID "<< nBXIds[iTrigClass] << ")"<< endl;
1284
1285         }
1286
1287       }
1288
1289       Float_t ratioToB[kNClasses];
1290       ratioToB[kClassE]  = nBXIds[kClassE]  >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassE]   : 0;
1291       ratioToB[kClassA]  = nBXIds[kClassA]  >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassA]   : 0;
1292       ratioToB[kClassC]  = nBXIds[kClassC]  >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassC]   : 0;
1293       ratioToB[kClassAC] = nBXIds[kClassAC] >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassAC]  : 0;
1294       Printf("Ratio between the BX ids in the different trigger classes:");
1295       Printf("  B/E  = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassE], ratioToB[kClassE] );
1296       Printf("  B/A  = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassA], ratioToB[kClassA] );
1297       Printf("  B/C  = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassC], ratioToB[kClassC] );
1298       Printf("  B/AC = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassAC],ratioToB[kClassAC]);
1299       Int_t nHistStat = 2;
1300
1301       // 1. loop over all cols
1302       for(Int_t iHistStat = 0; iHistStat < nHistStat; iHistStat++){
1303         Int_t ncol = fHistStatistics[iHistStat]->GetNbinsX();
1304         Float_t good1 = 0;
1305         for(Int_t icol = 1; icol <= ncol; icol++) {
1306           Int_t nEvents[kNClasses] = {0}; // number of events should be reset at every column
1307           // For all trigger classes, add up over row matching trigger mask (as selected before)
1308           for(Int_t iTrigClass = 0; iTrigClass < kNClasses; iTrigClass++){
1309             for(Int_t irow = 0; irow < nrows[iTrigClass]; irow++) {
1310               nEvents[iTrigClass] += (Int_t) fHistStatistics[iHistStat]->GetBinContent(icol,rows[iTrigClass][irow]);
1311             }
1312             //      cout << "Events " << classFlags[iTrigClass] << " ("<<icol<<") " << nEvents[iTrigClass] << endl;
1313           }
1314           if (nEvents[kClassB]>0) {
1315             Float_t acc  = ratioToB[kClassE]*nEvents[kClassE];
1316             Double_t acc_err = TMath::Sqrt(ratioToB[kClassE]*ratioToB[kClassE]*nEvents[kClassE]);
1317             //      Int_t bg   = cint1A + cint1C - 2*acc;
1318
1319              // Assuming that for a given class the triggers are either recorded as A+C or AC
1320             Float_t bg   = nEvents[kClassAC] > 0 ?
1321               fBIFactorAC*(ratioToB[kClassAC]*nEvents[kClassAC] - 2*acc):
1322               fBIFactorA* (ratioToB[kClassA]*nEvents[kClassA]-acc) +
1323               fBIFactorC* (ratioToB[kClassC]*nEvents[kClassC]-acc) ;
1324
1325             // cout << "-----------------------" << endl;
1326             // cout << "Factors: " << fBIFactorA << " " << fBIFactorC << " " << fBIFactorAC << endl;
1327             // cout << "Ratios: "  << ratioToB[kClassA] << " " << ratioToB[kClassC] << " " << ratioToB[kClassAC] << endl;
1328             // cout << "Evts:   "  << nEvents[kClassA] << " " << nEvents[kClassC] << " " << nEvents[kClassAC] << " " <<  nEvents[kClassB] << endl;
1329             // cout << "Acc: " << acc << endl;
1330             // cout << "BG: " << bg << endl;
1331             // cout  << "  " <<   fBIFactorA* (ratioToB[kClassA]*nEvents[kClassA]-acc) <<endl;
1332             // cout  << "  " <<   fBIFactorC* (ratioToB[kClassC]*nEvents[kClassC]-acc) <<endl;
1333             // cout  << "  " <<   fBIFactorAC*(ratioToB[kClassAC]*nEvents[kClassAC] - 2*acc) << endl;
1334             // cout << "-----------------------" << endl;
1335
1336             Float_t good = Float_t(nEvents[kClassB]) - bg - acc;
1337             if (icol ==1) good1 = good;
1338             //      Float_t errGood     = TMath::Sqrt(2*(nEvents[kClassA]+nEvents[kClassC]+nEvents[kClassE]));// Error on the number of goods assuming only bg fluctuates
1339             //      DeltaG^2 = B + FA^2 A + FC^2 C + Ratio^2 (FA+FC-1)^2 E.
1340             Float_t errGood     = nEvents[kClassAC] > 0 ?
1341               TMath::Sqrt( nEvents[kClassB] +
1342                            fBIFactorAC*fBIFactorAC*ratioToB[kClassAC]*ratioToB[kClassAC]*nEvents[kClassAC] +
1343                            ratioToB[kClassE] * ratioToB[kClassE] *
1344                            (fBIFactorAC - 1)*(fBIFactorAC - 1)*nEvents[kClassE]) :
1345               TMath::Sqrt( nEvents[kClassB] +
1346                            fBIFactorA*fBIFactorA*ratioToB[kClassA]*ratioToB[kClassA]*nEvents[kClassA] +
1347                            fBIFactorC*fBIFactorC*ratioToB[kClassC]*ratioToB[kClassC]*nEvents[kClassC] +
1348                            ratioToB[kClassE] * ratioToB[kClassE] *
1349                            (fBIFactorA + fBIFactorC - 1)*(fBIFactorA + fBIFactorC - 1)*nEvents[kClassE]);
1350
1351             Float_t errBG = nEvents[kClassAC] > 0 ?
1352               TMath::Sqrt(fBIFactorAC*fBIFactorAC*ratioToB[kClassAC]*ratioToB[kClassAC]*nEvents[kClassAC]+
1353                           4*ratioToB[kClassE]*ratioToB[kClassE]*(fBIFactorAC*fBIFactorAC)*nEvents[kClassE]) :
1354               TMath::Sqrt(fBIFactorA*fBIFactorA*ratioToB[kClassA]*ratioToB[kClassA]*nEvents[kClassA]+
1355                           fBIFactorC*fBIFactorC*ratioToB[kClassC]*ratioToB[kClassC]*nEvents[kClassC]+
1356                           ratioToB[kClassE]*ratioToB[kClassE]*(fBIFactorA+fBIFactorC)*(fBIFactorA+fBIFactorC)*nEvents[kClassE]);
1357
1358
1359             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowBG,bg);
1360             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowBG,errBG);
1361             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowAcc,acc);
1362             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowAcc,acc_err);
1363             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowGood,good);
1364             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowGood,errGood);
1365
1366 #ifdef VERBOSE_STAT
1367             //kStatRowBG=0,kStatRowAcc,kStatRowBGFrac,kStatRowAccFrac,kStatRowErrGoodFrac,kStatRowGoodFrac,kStatRowGood,kStatRowErrGood
1368             Float_t accFrac   = Float_t(acc) / nEvents[kClassB]  *100;
1369             Float_t errAccFrac= Float_t(acc_err) / nEvents[kClassB]  *100;
1370             Float_t bgFrac    = Float_t(bg)  / nEvents[kClassB]  *100;
1371             Float_t goodFrac  = Float_t(good)  / good1 *100;
1372             Float_t errGoodFrac = errGood/good1 * 100;
1373             Float_t errFracBG = bg > 0 ? TMath::Sqrt((errBG/bg)*(errBG/bg) + 1/nEvents[kClassB])*bgFrac : 0;
1374             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowBGFrac,bgFrac);
1375             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowBGFrac,errFracBG);
1376             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowAccFrac,accFrac);
1377             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowAccFrac,errAccFrac);
1378             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowGoodFrac,goodFrac);
1379             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowErrGoodFrac,errGoodFrac);
1380             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowErrGood,errGood);
1381 #endif
1382           }
1383         }
1384       }
1385       for (Int_t iTrigClass = 0; iTrigClass < kNClasses; iTrigClass++){
1386         delete [] rows[iTrigClass];
1387       }
1388     }
1389   }
1390
1391   fHistStatistics[0]->Write();
1392   fHistStatistics[1]->Write();
1393   if(fHistBunchCrossing ) fHistBunchCrossing ->Write();
1394   if(fHistTriggerPattern) fHistTriggerPattern->Write();
1395
1396   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
1397   for (Int_t i=0; i < count; i++)
1398   {
1399     TString triggerClass = "trigger_histograms_";
1400     if (i < fCollTrigClasses.GetEntries())
1401       triggerClass += ((TObjString*) fCollTrigClasses.At(i))->String();
1402     else
1403       triggerClass += ((TObjString*) fBGTrigClasses.At(i - fCollTrigClasses.GetEntries()))->String();
1404
1405     gDirectory->mkdir(triggerClass);
1406     gDirectory->cd(triggerClass);
1407
1408     static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i))->SaveHistograms();
1409
1410     gDirectory->cd("..");
1411   }
1412
1413   if (fBackgroundIdentification)
1414   {
1415     gDirectory->mkdir("background_identification");
1416     gDirectory->cd("background_identification");
1417
1418     fBackgroundIdentification->GetOutput()->Write();
1419
1420     gDirectory->cd("..");
1421   }
1422
1423   if (folder)
1424     gDirectory->cd("..");
1425 }
1426
1427 Int_t AliPhysicsSelection::GetStatRow(const char * triggerBXClass, UInt_t offlineTriggerType, UInt_t ** rowIDs) const {
1428   // Puts inside the array rowIDs the row number for a given offline
1429   // trigger in a given bx class. Returns the total number of lines
1430   // matching the selection
1431   // triggerBXClass can be either "A", "AC", "B" or "E"
1432   // offlineTriggerType is one of the types defined in AliVEvent
1433   // User should delete rowIDs if no longer needed
1434
1435   if(!fHistStatistics[0]) {
1436     AliWarning("Not initialized, returning 0");
1437     return 0;
1438   }
1439   const Int_t nrows = fHistStatistics[0]->GetNbinsY();
1440
1441   // allocate memory for at maximum nrows
1442   Int_t nMatches = 0;
1443   (*rowIDs) = new UInt_t[nrows];
1444
1445   // Build regular expression. look for a +, followed by the beginning
1446   // of a word. Within the word, look for the class id after any
1447   // number of any char, but at most one dash ("[^-]*-?"), followed by
1448   // a - and then any char (".*") and at the class id ("\\&(\\d)")
1449   // The class id is stored.
1450   // WARNING: please check this if the trigger classes change
1451   TPRegexp re(Form("\\+\\b[^-]*-?%s-.*\\&(\\d)",triggerBXClass));
1452   // Loop over rows and find matching ones:
1453   for(Int_t irow = 1; irow <= nrows; irow++){
1454     TObjArray * matches = re.MatchS(fHistStatistics[0]->GetYaxis()->GetBinLabel(irow));
1455     if (matches->GetEntries()) {
1456       TString s = ((TObjString*)matches->At(1))->GetString();
1457       if(UInt_t(s.Atoi()) & offlineTriggerType) { // bitwise comparison with the requested mask
1458         //      cout << "Marching " << s.Data() << " " << offlineTriggerType << " " << fHistStatistics[0]->GetYaxis()->GetBinLabel(irow) << endl;
1459         (*rowIDs)[nMatches] = irow;
1460         nMatches++;
1461       }
1462     }
1463     delete matches;
1464   }
1465
1466   return nMatches;
1467 }
1468
1469 void AliPhysicsSelection::SetBIFactors(const AliESDEvent * aESD) {
1470   // Set factors for realtive bunch intesities
1471   if(!aESD) {
1472     AliFatal("ESD not given");
1473   }
1474   Int_t run = aESD->GetRunNumber();
1475   if (run > 105268) {
1476     // intensities stored in the ESDs
1477     const AliESDRun* esdRun = aESD->GetESDRun();
1478     Double_t intAB = esdRun->GetMeanIntensityIntecting(0);
1479     Double_t intCB = esdRun->GetMeanIntensityIntecting(1);
1480     Double_t intAA = esdRun->GetMeanIntensityNonIntecting(0);
1481     Double_t intCC = esdRun->GetMeanIntensityNonIntecting(1);
1482
1483     // cout << "INT " <<intAB <<endl;
1484     // cout << "INT " <<intCB <<endl;
1485     // cout << "INT " <<intAA <<endl;
1486     // cout << "INT " <<intCC <<endl;
1487
1488     if (intAB > -1 && intAA > -1) {
1489       fBIFactorA = intAB/intAA;
1490     } else {
1491       AliWarning("Cannot set fBIFactorA, assuming 1");
1492     }
1493
1494     if (intCB > -1 && intCC > -1) {
1495       fBIFactorC = intCB/intCC;
1496     } else {
1497       AliWarning("Cannot set fBIFactorC, assuming 1");
1498     }
1499
1500     if (intAB > -1 && intAA > -1 &&
1501         intCB > -1 && intCC > -1) {
1502       fBIFactorAC = (intAB+intCB)/(intAA+intCC);
1503     } else {
1504       AliWarning("Cannot set fBIFactorAC, assuming 1");
1505     }
1506
1507   }
1508   else {
1509     // First runs. Intensities hardcoded
1510     switch(run) {
1511     case 104155:
1512       fBIFactorA = 0.961912722908;
1513       fBIFactorC = 1.04992336081;
1514       break;
1515     case 104157:
1516       fBIFactorA = 0.947312854998;
1517       fBIFactorC = 1.01599706417;
1518       break;
1519     case 104159:
1520       fBIFactorA = 0.93659320151;
1521       fBIFactorC = 0.98580804207;
1522       break;
1523     case 104160:
1524       fBIFactorA = 0.929664189926;
1525       fBIFactorC = 0.963467679851;
1526       break;
1527     case 104315:
1528       fBIFactorA = 1.08939104979;
1529       fBIFactorC = 0.931113921925;
1530       break;
1531     case 104316:
1532       fBIFactorA = 1.08351880974;
1533       fBIFactorC = 0.916068345845;
1534       break;
1535     case 104320:
1536       fBIFactorA = 1.07669281245;
1537       fBIFactorC = 0.876818744763;
1538       break;
1539     case 104321:
1540       fBIFactorA = 1.00971079602;
1541       fBIFactorC = 0.773781299076;
1542       break;
1543     case 104792:
1544       fBIFactorA = 0.787215863962;
1545       fBIFactorC = 0.778253173071;
1546       break;
1547     case 104793:
1548       fBIFactorA = 0.692211363661;
1549       fBIFactorC = 0.733152456667;
1550       break;
1551     case 104799:
1552       fBIFactorA = 1.04027825161;
1553       fBIFactorC = 1.00530825942;
1554       break;
1555     case 104800:
1556       fBIFactorA = 1.05309910671;
1557       fBIFactorC = 1.00376801855;
1558       break;
1559     case 104801:
1560       fBIFactorA = 1.0531231922;
1561       fBIFactorC = 0.992439666758;
1562       break;
1563     case 104802:
1564       fBIFactorA = 1.04191478134;
1565       fBIFactorC = 0.979368585208;
1566       break;
1567     case 104803:
1568       fBIFactorA = 1.03121314094;
1569       fBIFactorC = 0.973379962609;
1570       break;
1571     case 104824:
1572       fBIFactorA = 0.969945926722;
1573       fBIFactorC = 0.39549745806;
1574       break;
1575     case 104825:
1576       fBIFactorA = 0.968627213937;
1577       fBIFactorC = 0.310100412205;
1578       break;
1579     case 104841:
1580       fBIFactorA = 0.991601393212;
1581       fBIFactorC = 0.83762204722;
1582       break;
1583     case 104845:
1584       fBIFactorA = 0.98040863886;
1585       fBIFactorC = 0.694824205793;
1586       break;
1587     case 104867:
1588       fBIFactorA = 1.10646173412;
1589       fBIFactorC = 0.841407246916;
1590       break;
1591     case 104876:
1592       fBIFactorA = 1.12063452421;
1593       fBIFactorC = 0.78726542895;
1594       break;
1595     case 104890:
1596       fBIFactorA = 1.02346137453;
1597       fBIFactorC = 1.03355663595;
1598       break;
1599     case 104892:
1600       fBIFactorA = 1.05406025913;
1601       fBIFactorC = 1.00029166135;
1602       break;
1603     case 105143:
1604       fBIFactorA = 0.947343384349;
1605       fBIFactorC = 0.972637444408;
1606       break;
1607     case 105160:
1608       fBIFactorA = 0.908854622177;
1609       fBIFactorC = 0.958851103977;
1610       break;
1611     case 105256:
1612       fBIFactorA = 0.810076150206;
1613       fBIFactorC = 0.884663561883;
1614       break;
1615     case 105257:
1616       fBIFactorA = 0.80974912303;
1617       fBIFactorC = 0.878859123479;
1618       break;
1619     case 105268:
1620       fBIFactorA = 0.809052110679;
1621       fBIFactorC = 0.87233890989;
1622       break;
1623     default:
1624       fBIFactorA = 1;
1625       fBIFactorC = 1;
1626     }
1627   }
1628
1629 }