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